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Summary 


Turbulent wall-bounded complex flows are commonly encountered in engineering practice 
and are of considerable interests in a variety of industrial applications. The presence of a wall 
significantly affects turbulence characteristics. In addition to the wall effects, turbulent wall- 
bounded flows become more complicated by the presence of additional body forces (e.g. 
centrifugal force and Coriolis force) and complex geometry. Most near-wall Reynolds stress 
models are developed from a high-Reynolds-number model which assumes turbulence is 
homogenous (or quasi-homogenous). Near-wall modifications are proposed to include wall 
effects in near-wall regions. In this process, wall normals are introduced. Good predictions could 
be obtained by Reynolds stress models with wall normals. However, ambiguity arises when the 
models are applied in flows with multiple walls. 

Many models have been proposed to model turbulent flows. Among them, Reynolds stress 
models, in which turbulent stresses are obtained by solving the Reynolds stress transport 
equations, have been proved to be the most successful ones. To apply the Reynolds stress models 
to wall-bounded flows, near-wall corrections accounting for the wall effects are needed, and the 
resulting models are called near-wall Reynolds stress models. In most of the existing near-wall 
models, the near-wall corrections invoke wall normals. These wall-dependent near-wall models 
are difficult to implement for turbulent flows with complex geometry and may give inaccurate 
predictions due to the ambiguity of wall normals at comers connecting multiple walls. 

The objective of this study is to develop a more general and flexible near-wall Reynolds 
stress model without using any wall-dependent variable for wall-bounded turbulent flows. With 
the aid of near-wall asymptotic analysis and results of direct numerical simulation, a new near-wall 
Reynolds stress model (NNWRS) is formulated based on Speziale et al.'s high-Reynolds-stress 
model with wall-independent near-wall corrections. Moreover, only one damping function is used 
for flows with a wide range of Reynolds numbers to ensure that the near-wall modifications 
diminish away from the walls. 


vii 



Various simple and complex turbulent wall-bounded flows are used to validate the NNWRS 
model. Model predictions agree reasonably well with available data from experiments, direct 
numerical simulation, or large eddy simulation. Complex flow features caused by the centrifugal 
force and Coriolis force as in swirling pipe flow, axially rotating pipe flow and channel flow with 
spanwise rotation are essentially captured by the model. The model is able to reproduce 
complicated flow phenomena induced by complex geometry, such as flow recirculation, 
reattachment and boundary-layer redevelopment in backward-facing step flow and secondary flow 
in three-dimensional square duct flow. In simple flows, including fully developed channel/pipe 
flow, Couette flow and boundary-layer flow, the wall effects are dominant, and the NNWRS 
model predicts less degree of turbulent anisotropy in the near-wall region compared with the wall- 
dependent near- wall Reynolds stress model (NWRS) developed by So et al.. However, in 
complex flows where other effects become more important, the NNWRS model performs almost 
as well as the NWRS model, and the advantages of the wall-independent model become more 
evident. The comparison of the predictions given by the two models also rectifies the 
misconception that the overshooting of skin friction coefficient in backward-facing step flow 
prevalent in those near-wall models with wall normals is caused by the use of wall normals. 
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1 INTRODUCTION 


1 . 1 Background 

Turbulence closure modeling to solve the mean-averaged Navier-Stokes equations has been 
pursued for nearly a century. In the last several decades, with the great advances in computer 
science, numerous turbulence models have been developed to simulate and predict more and more 
complex turbulent flows. In most engineering problems, flow fields are bounded by walls. The 
presence of a wall affects significantly the turbulence characteristics: it reduces the length scale of 
the momentum fluctuations and increases the dissipation rate; it reflects the pressure fluctuation, 
thereby inhibiting the transfer of turbulence energy into fluctuations normal to the wall; it enforces 
a no-slip boundary condition, thus ensuring that within a wall-adjacent sublayer, however thin, 
turbulent stress is negligible and viscous effect on transport processes becomes vitally important 
(Launder 1989). 

In general, there are two types of turbulence models: eddy-viscosity models and second- 
order closure models (Reynolds stress models). In the eddy-viscosity models, the Reynolds stress 
tensor is assumed to be the product of an eddy viscosity and the mean strain-rate tensor. In the 
second-order closure models, turbulent Reynolds stresses are obtained by solving the Reynolds 
stress transport equations and no relation is presumed between the Reynolds stresses and the mean 
strain rates. As a result, the second-order closure models are more general and can be used to 
simulate turbulence anisotropy and the turbulent stress redistribution process in the near-wall 
region. Furthermore, other important factors such as streamline curvature, centrifugal and Coriolis 
forces could be properly accounted for in the Reynolds stress models, whereas they are rather 
difficult to be included in the eddy-viscosity models. 

Practically all turbulence closure models invoke the large Reynolds number assumption, thus 
allowing the viscous effect to be neglected as a first approximation. Consequently, the high- 
Reynolds-number models cannot be applied to the near-wall region where viscous effect is 
dominant. To remedy this drawback, wall-functions have been used to bridge the gap between the 
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wall and the location beyond which the high-Reynolds-number models are applicable. The wall- 
functions are derived under the assumptions of equilibrium turbulence and constant shear stress 
near a wall, which are less likely to be valid for complex turbulent shear flows. Therefore, many 
attempts have been made to develop near-wall second-order closure models. 

As a wall is approached, the local Reynolds number decreases and the intensity of the 
anisotropy of the turbulence field increases due to the reflection of the fluctuating pressure by the 
wall (Kim 1989). This effect is commonly known as wall blocking (wall-reflection or pressure- 
echo) effects. Although the influence of the viscous effect and wall blocking effect on the turbulent 
motion are very different by nature, these effects have been frequently modeled together in near- 
wall turbulence modeling. All the existing near-wall Reynolds-stress models are developed based 
on the high-Reynolds-number models. To account for wall effects, in principle, near-wall 
modifications for the turbulent diffusion, velocity-pressure-gradient correlation, and dissipation 
rate tensor are needed so that these terms have the correct asymptotic behavior near the wall. 
Asymptotic analysis shows that the turbulent diffusion term is three orders of magnitude smaller 
than the velocity-pressure-gradient correlation and the dissipation rate tensor in the near-wall 
region. Therefore, no near-wall modification is actually required for the turbulent diffusion term. 
On the other hand, the velocity-pressure-gradient correlation and dissipation rate tensor are order 
one quantities in the near-wall region, and near-wall modifications for these terms are necessary 
for the models to be applicable in the near-wall region. To ensure that the near-wall modifications 
diminish away from the wall, a damping function analogous to that proposed by Van Driest (1956) 
is introduced so that the high-Reynolds-number models are recovered in the region far away from 
the wall. 

The velocity-pressure-gradient correlation and dissipation rate tensor play crucial roles in 
near-wall turbulence and represent different physical processes. The velocity-pressure-gradient 
correlation is usually partitioned into a pressure-strain part and a pressure diffusion part. Few 
models have explicitly considered the pressure diffusion effect (e.g. Lumley 1978). 
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Pressure diffusion is assumed to be either included into the modeling of the turbulent diffusion 
term, or insignificant in near-wall turbulence as it is treated when the Reynolds number is large. 
On the other hand, the contribution of the pressure-strain correlation is significant because it 
redistributes energy among different turbulent stress components. The pressure-strain correlation 
is a linear function of the pressure fluctuation, and according to the structure of the solution to the 
Poisson equation for the pressure fluctuation, it consists of a 'return' term, a 'rapid' term and a 
Stokes term representing the pressure reflection caused by the wall (Mansour et al. 1988). Due to 
the difficulty to include the Stokes part in the model if local homogeneous turbulence is assumed, 
the Stokes term is neglected in most of the high-Reynolds-number models (e.g. Launder et al. 
1975). When these models are extended to the near-wall region, an additional term proposed to 
simulate the Stokes term is added back to the pressure-strain correlation (Prud'homme and 
Elghobashi 1983; Kebede et al. 1984; Shima 1988; Shih and Mansour 1990; Lai and So 1990). 
Another approach is to modify the coefficients of the return and rapid terms to justify the neglect of 
the Stokes term (Launder and Shima 1989). In some models, the Stokes term is absorbed into 
either the return part or the rapid part (Speziale et al. 1991; Launder and Tselepidakis 1993; 
Launder and Li 1994). To ensure that the near-wall corrections decay exponentially away from 
the wall, a function depending on the wall normal distance is usually introduced, and the resulting 
models are no longer coordinate-frame invariant. 

The dissipation rate tensor £„• in the Reynolds-stress equation is formulated as the 
combination of an isotropic solenoidal part and an anisotropic part according to Hanjalic and 
Launder (1976). The isotropic part is obtained by solving a transport equation for the dissipation 
rate £ of the kinetic energy. Hanjalic and Launder (1976) proposed the first £ -equation 
compatible with their Reynolds-stress model. In spite of its oversimplified form compared with 
the exact governing equation for dissipation rate, it performs reasonably well and becomes a 
prototype for later modifications. Modifications have been made mainly to reproduce the near- 
wall behavior of dissipation rate predicted by direct numerical simulation: a maximum value of £ 
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is expected to occur at the wall. To reproduce such behavior, near-wall corrections are introduced 
to the dissipation rate equation and wall dependent variables (wall normals) are used in damping 
functions (Shima 1988) or transformed dissipation terms (Lai and So 1990). The anisotropic 
correction part is constructed to satisfy the following two requirements: (i) the contraction of the 
dissipation rate tensor should be equal to two times the dissipation rate of the kinetic energy, i.e. 
££;, = 2e; (ii) the components of the dissipation rate tensor should satisfy the kinematic 
constraints proposed by Launder and Reynolds (1983). Several models have been proposed to 
improve the dissipation rate tensor in the near-wall region (e.g. Prud'homme and Elghobashi 1983; 
Kebede et al. 1984; So and Yoo 1986; Lai and So 1990). Among them, Lai and So's (1990) 
model satisfies the kinematic constraints and contraction requirement with the introduction of the 
wall unit vectors. Most recently, Shima (1995) proposed a three-term dissipation rate model. This 
model also satisfies both requirements but eliminates all wall normal unit vectors. 

In view of the above discussion, it can be concluded that most of the existing near-wall 
Reynolds-stress models explicitly invoke wall-dependent variables such as wall normal unit vector 
and wall normal distance to account for the viscous and wall effects. The influence of each wall in 
this wall-dependent model is also assumed to act independently. For turbulence flows with 
complex geometry, which are commonly encountered in most engineering problems, the use of 
wall normals at corners connecting multiple walls becomes ambiguous. Therefore, in this 
situation, these models are difficult to implement and may introduce inaccuracy in the prediction of 
complex turbulent flow field. The present study is aimed to improve the existing near-wall 
Reynolds-stress models for complex flows caused by complex geometry, as well as streamline 
curvature and additional body forces such as centrifugal and Coriolis forces in flows with swirl 
and rotation. 

1 . 2 Present Objectives 

The present research is to develop a near-wall Reynolds-stress model for complex turbulent 
flows without introducing any wall-dependent variables (wall normals). Several near-wall second- 
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order closure models have been proposed at Arizona State University (ASU) in the last decade. 
They are So and Yoo (1987), Lai and So (1990), So et al. (1991) and So et al. (1994). With the 
aid of asymptotic analysis, the near-wall modifications in these models were developed based on 
the high-Reynolds-number model of Launder et al. (1975) except that of So et al's (1994), which 
was based on the model of Speziale et al. (1991). This latest model (hereafter we will call it 
NWRS model) has been validated against various applications ranging from incompressible fully 
developed channel flow to high Mach number compressible boundary-layer flow in a wide range 
of Reynolds number. Although these models developed by the ASU group are asymptotically 
consistent to a certain extent in the near-wall region, all of them are dependent on wall variables in 
one way or another. The present study is part of the on-going near-wall second-order closure 
model development project. An attempt to develop a near-wall Reynolds-stress model without wall 
dependence has been carried out by Launder and Li (1994) with partial success. The wall normals 
were removed from the pressure-strain correlation in their Reynolds-stress model, but remain in 
the dissipation rate tensor. Therefore, the model developed in this study is expected to be the first 
near-wall Reynolds-stress model with no wall dependence. The challenge of near-wall modeling 
lies on the lack of thorough understanding of wall effects on turbulence. A near-wall second-order 
closure model without wall-dependent variables poses more difficulties because of the limited 
number of solvable variables and suitable formulations. Guided by the results of direct numerical 
simulation, we will develop an asymptotically correct near-wall model without wall normals in this 
study. The new model is a modification of NWRS which can accommodate flow complexities 
such as streamline curvature and additional body forces related to fluid rotation. The new model, 
designated as NNWRS, is expected to perform as well as the wall-dependent NWRS model in a 
wide range of complex flow conditions. 

1 . 3 Report Outline 

In Section 2, a new near-wall Reynolds-stress (NNWRS) model without wall normals is 
developed based on the pressure-strain model of Speziale et al. (1991). Asymptotic analysis and 
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results from direct numerical simulation are used to guide the near-wall modifications for the 
pressure-strain correlation, dissipation rate tensor, and dissipation rate equation. For comparison, 
the formulation for the near-wall Reynolds stress (NWRS) model with wall normals proposed by 
So et al. (1994) is also presented, and its predictions will be compared with those given by the 
NNWRS model in Sections 3, 4 and 5, together with available data from experiments, direct 
numerical simulation and large eddy simulation. 

In Section 3, simple internal and external flows with a wide range of Reynolds numbers are 
first used to validate the new proposed model. These flows include fully-developed channel/pipe 
flow, Couette flow, and boundary-layer flow with zero pressure gradient. In addition to assessing 
its performance on the mean and turbulence fields, the model's ability to replicate Reynolds 
number effects on the mean flow and second-order statistics is also examined. 

In Section 4, the NNWRS model is applied to swirling pipe flow, axially rotating pipe flow, 
and a channel flow with spanwise rotation to demonstrate its ability to reproduce complicated flow 
phenomena caused by streamline curvature, centrifugal force, and Coriolis force in a wide range of 
Reynolds numbers, swirl numbers, and rotation numbers. 

In Section 5, two flows with multiple walls, a two-dimensional backward-facing step flow 
and a three-dimensional square duct flow, are used to test the model’s ability to replicate flow 
phenomena resulting from complex geometry, such as flow recirculation, reattachment, and 
boundary-layer redevelopment in backward-facing step flow, and the secondary flow induced by 
turbulence in square duct flow. The advantages of the NNWRS model for flows with complex 
geometry are demonstrated. 

Section 6 summarizes the performance of the proposed near-wall Reynolds-stress model 
without wall normals and presents the conclusions drawn. 
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2. NEAR-WALL REYNOLDS-STRESS MODELING 


2.1 Mean Flow Equations 

The present study considers turbulent flows of a viscous, incompressible fluid. The fluid 
motion with the presence of rotation is described by the Navier-Stokes equations, which can be 
written in Cartesian form or notation as: 



1 dp d 2 u i _ _ _ 

~p^ +V ~d^~ eijkQjUk 


where w, is the instantaneous velocity vector, p is the modified pressure including the centrifugal 
force potential, Qj is the angular velocity, e ljk is the permutation tensor, and p and v are the fluid 

density and viscosity, respectively. Equations (2-1) and (2-2) express conservation of mass and 
momentum per unit mass, respectively. 

The N-S equations (1) and (2) are fundamental governing equations and can be applied to 
both laminar and turbulent flows. With specified boundary and initial conditions, in principle, (1) 
and (2) can be solved numerically. However, for turbulent flows, the resolution of the small scale 
turbulent fluctuations requires very fine grids; as a result, direct numerical simulations (DNS) of 
turbulent flows at high Reynolds number are extremely difficult, if not impossible. An alternative 
approach is to study the mean flow field with the consideration of the influence of turbulence. The 
mean flow field can be obtained through time, spatial or ensemble averaging. Time averaging has 
been widely used in engineering since stationary turbulent flows are most frequently encountered. 
This method is also adopted in this study. By time averaging, the velocity and pressure fields are 
decomposed into a mean and a fluctuating part: 

m, = U t + «, , p= P + p , (2-3) 


where capital letters denote the time-averaged mean quantities and small case symbols represent 
the fluctuating quantities. 
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Substituting (2-3) into (2-1) and (2-2) and taking the time averaging of the resulting 
equations, we obtain the governing equations for the mean flow field, or 


dU 1 dP d 2 U Bum. 

■ <«> 

Note that the effect of turbulence on the mean flow is represented by the Reynolds stress term u,Uj 

(i.e. correlations between different fluctuating velocity components) on the right hand side of (2- 
5). Equations (2-4) and (2-5) are called the Reynolds equations for the mean flow field. 

To solve the Reynolds equations, one must relate the Reynolds stress term to the mean 
velocity to close the equations. In the following section we derive the transport equations for the 
Reynolds stresses and for the dissipation rate. 

2.2 Time-Averaged Turbulence Transport Equations 

The transport equations for the Reynolds stresses and the dissipation rate can be derived from 
equations (2-1) and (2-2). Defining a Navier-Stokes operator 

.... . du. . du, 1 dp d 2 U: _ „ _ ._ 

Nlu ‘ )= iF +u ^ + p&r v ^ ,Uk • ( ’ 


and carrying out the following time averaging 
u j N(u j ) + UjNi&j ) = 0 


we obtain the Reynolds-stress transport equations: 

diijiij 3 ( diijUj \ dujUjU^ f dU 


Bxl Bxl 1 dx 


+ -«i u k 




. _ / \ f Bu. duj 

— 2 Gif, yi jU in i > l i :m + U l li nl € jj <)n J + p ^ 
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The left-hand side is the convection term and the terms on the right-hand side of (2-8) are viscous 
diffusion, turbulent diffusion, production due to mean velocity strain, production due to rotation, 
pressure-strain correlation, pressure diffusion and viscous dissipation. Similarly, the equation for 
the dissipation rate e, defined as 


e = 


du. du, 

V d^d^ ’ 


(2-9a) 


can be derived from 


2 v^- ‘ - 3 [AT(m,)] = 0 , 


dxj d.Xj 


(2-9b) 


where JV(w ( )is the Navier-Stokes operator defined in (2-6) with w, p p. After tedious 
algebraic manipulation, the following equation for e is obtained 


U, 


de 


dxi dx l 


V 


de 


dx 


+ - 


kj 


dx L 


-VUi, 


du , du ; 


dx m dx m j 


2v 


dui duj dU t 
dxu dxu dx. 


_ du k du,, dU : „ 

— 2v — — . L -2v\ 


-2 


V d 


dx i dxj dxj 

r 


du l 

dx 


J J 


2 v du t du t du k 


dx k dx: 


dx k dx m dx„ 


p dxj 


dp duj 
dx m dx m j 


- 2 v 


2 d 2 U. 


OU: 


dx k dx m dx k dx m 


( 2 - 10 ) 


The terms on the right-hand side of (2-10) are generally regarded as viscous diffusion, turbulent 
diffusion, mixed production, production by mean velocity gradient, gradient production, 
turbulence production, pressure diffusion and turbulent destruction of the dissipation. 

To facilitate the discussions on each term in the Reynolds-stress and dissipation rate transport 
equations, we recast (2-8) and (2-10) in the following symbolic forms: 

Qj - dH + d? + p i j + R ij +n ’ J - , (Mi) 

C e = D v t +Dl + P' c + li + Pl + P* + Dr-y , ( 2 - 12 ) 

respectively, where each symbol represents each position-corresponding term on the right-hand 
side of the respective full equations. Note that in (2-11) the turbulent diffusion term is 
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(2- 13a) 



(2- 13b) 

(2- 13c) 
(2- 13d) 

(2-13e) 

(2-14) 


need modeling if the Reynolds equations for the mean flow field are to be closed. 

2.3 High-Reynolds-Number Modeling 

In this section, we discuss how the turbulent diffusion Dj, velocity-pressure-gradient 
correlation n* and dissipation rate tensor in (2-11) are modeled under the assumption of high 

Reynolds number (the resulting models are called the high-Reynolds-number models). For high 
Reynolds number flows, the viscous effect can be neglected. Irrespective of the analytical 
arguments employed for modeling, all known high-Reynolds-number models use the Reynolds 
stress gradients to express the turbulent diffusion term (some models also employ the gradient of 
the turbulence scale or the gradient of the scale-supplying variable, such as the dissipation rate £). 
Five formulations have been proposed for the turbulent diffusion term and they are given by Daly 
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and Harlow (1970), Shir (1973), Hanjalic and Launder (1972), Mellor and Herring (1973) and 
Lumley (1978). The invariant form of turbulent diffusion by Hanjalic and Launder (1972) 

d 


dx L 


C - 
S £ 


u,u 


i u l ■ 


V 


dujU k 

dx> 


+ U;U 


J u l ■ 


du k Uj 

dx. 


+ u k u l 


dn t u ■ x 


dx 


l ) 


(2-15) 


is adopted for the present model development, where k = 1 / 2 u i u i is the turbulence kinetic energy 
and C s is a constant. This formula is not only tensorily consistent with the exact expression but 

also found to perform better in several types of turbulent flows. 

Velocity-pressure-gradient correlation n* is traditionally partitioned into a pressure-strain term 
n (> and a pressure diffusion termD? (2- 13b). The contribution from is usually neglected for 

high-Reynolds-number flows or is argued to be included in turbulent diffusion model (2-15). In 
either case, the velocity-pressure-gradient term n* is considered to be the same as the pressure- 
strain term Tly in high-Reynolds-number modeling. 


For incompressible flows, the pressure fluctuation is governed by the following Poisson 
equation 


P 



dxj dx t 


<9 2 

dxjdxj 


( 


UjU j - Ujllj 



(2-16) 


with boundary condition 

dp _ d 2 v 
dy dy 2 


(2-17) 


where v is the velocity fluctuation in the wall normal direction (y direction). The Poisson 
equation and its boundary condition (2-17) are linear with respect to the fluctuating pressure p . 
Therefore, its solution can be splitted into three parts, a 'return' part, a 'rapid' part and a 'Stokes' 
part (Mansour et al. 1988). The return pressure, p ^ , is defined as the solution of the following 
problem 


P 



<9 2 

dx t dxj 


[uiUj-UjUj] 


(2- 18a) 


with the boundary condition at the wall given by 
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3a. = o 

dy 


(2- 18b) 


The rapid pressure, pj, is defined as the solution to 



dU, duj 
dxj dx i 


(2- 19a) 


with the boundary condition at the wall given by 

-^1 = 0 (2- 19b) 

ft 

Finally, the Stokes pressure, p s , is defined as the solution to 

-V 2 Ps = 0 , (2-20a) 

P 

with the boundary condition at the walls specified as 

= (2-20b) 

By dy 2 

This split resolves the question of whether to add the inhomogeneous boundary condition to the 
return part of the pressure or to the rapid part. It does not affect the wall effect on the pressure 
fluctuations. The pressure-strain term is linear in p (2- 13c) and accordingly the Stokes pressure- 
strain statistics can be added to either the rapid pressure-strain term or the return pressure-strain 
term without affecting the wall effect on the pressure-strain correlation. However, the Stokes 
pressure-strain term (representing the wall effect) is neglected in most models due to the 
mathematical difficulty and the lack of understanding of the wall effect on turbulent flow field. 
Launder et al.'s (1975) pressure-strain correlation (hereafter denoted as the LRR model) has been 
widely accepted for most Reynolds-stress models because of its good performance. In the LRR 
model, a wall reflection term was introduced to simulate the Stokes pressure-strain term and to 
compensate for the stress anisotropy due to the presence of the wall. The inclusion of the wall 
reflection term has been regarded as indispensable to simulate the wall effects. However, it causes 
a major drawback since this term involves the distance from the wall and is not coordinate frame 
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invariant. On the other hand, the SSG model (Speziale et al. 1991) incorporated the Stokes term 
into the return and rapid parts, and no additional wall reflection term is proposed. 

The recent trend for pressure-strain correlation modeling is to include additional nonlinear 
terms to the conventional linear terms. The inclusion of nonlinear terms brings more flexibility to 
satisfy certain kinematic constraints (e.g. realizability). Among the nonlinear models, the SSG 
model retains the linear parts in the LRR model, but with some coefficients depending on the 
turbulent stress invariants and turbulence production (hence, it is a quasi-linear model). Other 
models contain the quadratic and cubic terms in the rapid part and therefore are much more 
complicated to use in engineering applications (Choi and Lumley 1984; Shih and Lumley 1985; 
Craft and Launder 1991; Launder and Tselepidikis 1993). 

With appropriate rearrangement and transformation (Appendix A), the SSG model can be 
rewritten in the following form: 

n y =-(Qe + C,‘ P)b i} + c 2 e(b ik b kj - ^n<5 (y ) + C 5 k(b jm e lkm + b im e jkm )L\ 

-a i (P ij ^PS i j)-p ] (D ij -^PS i j)-2(y l +^n m )kS ij , (2-21) 

where C] Cj\ C 2 , C 3 , C 5 , cq, and are constant coefficients, and the expressions for 

P, b t j, FI, Pjj, Dy and are given in Appendix A. By rewriting the SSG model in this form, the 

£ $ 

LRR model is readily recovered by setting Q = C 2 = C 3 = C 5 =0, and the meaning of each 

individual term in (2-21) can be easily identified. The first two terms on the right-hand side of (2- 
21) (Cj* andC 2 terms) are the nonlinear return part, the third term (C 5 term) is to account for 

rotation effect, and the rest ( cq, /)| and y l terms) are similar to the rapid part in the LRR model. 
The SSG model gives better predictions than the LRR model does for flows with streamline 
curvature and rotation (Speziale et al. 1992). In view of the above discussion, the SSG model is 
selected as the base for the development of near-wall second-order closure models in this study. 

In high-Reynolds-number turbulence, the dissipation rate tensor is assumed to be isotropic 
and the form proposed by Kolmogorov ( 1 94 1 ) is often adopted, which is 
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( 2 - 22 ) 



where the isotropic dissipation rate £ is obtained by solving the following modeled transport 
equation. 


U t 


de 


dxi dxi 


de 

V dx k J 


+ • 


dxi. 


C r — UjU, 


i u k 


de_ 

dx: 


i J 


+ C el f t ^p-C e2 f 2 j 


(2-23) 


and C E C £l and C e2 are constant coefficients, /j and f 2 are the damping functions. 


2.4 Near- Wall Reynolds-Stress Modeling Based on the SSG Model 

The presence of a wall affects significantly the turbulence behaviors in many aspects. Near- 
wall turbulence modeling is to simulate the wall effects (including viscous and blocking effects) on 
turbulent flows. In this section, the high-Reynolds-number models discussed in the previous 
section are extended to the near-wall region by incorporating asymptotically correct near-wall terms 
to the model of the pressure-strain correlation, the dissipation rate transport equation and the model 
of the dissipation rate tensor. The derived near-wall modifications are ensured to vanish away 
from the wall through the use of damping functions. Two near-wall Reynolds-stress models will 
be presented in this section, one is wall-dependent (hereafter referred to as NWRS), the other is 
wall-independent (hereafter referred to as New NWRS or NNWRS). The wall-dependent model 
has been reported in So et al.(1994), in which the author is one of the co-authors. 

Following Lai and So (1991), we expand the fluctuating velocity components, which satisfy 
the no-slip boundary conditions at the wall, into Taylor series in the near-wall region as 


2 3 

u = a]y + a 2 y + a 3 y + .... , 

(2-24a) 

v = b^y + b 2 y 2 + b^y 3 + .... , 

(2-24b) 

2 3 

w = qy + c 2 y +c 3 y + .... , 

(2-24c) 


where = (w,v,w), x, = (x,y,z), and the y-axis is taken to be normal to the wall, the x-axis is in 
the stream direction and the z-axis is normal to the (x, y)-plane. The coefficients a h b h c i are 
random functions of time, x, and z, but not y. For incompressible flows, b : = 0 is required to 
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satisfy the continuity equation. The near-wall behavior of each term in the Reynolds stress 
equation (2-11) can therefore be analyzed by the substitution of (2-24), except the behavior of the 
velocity-pressure-gradient term, which can be obtained by rearranging ( 2 - 11 ) as 
n * = C,y - Djj - D]j - Pjj - R { j + £,y. Table 2.1 shows the near-wall behavior of each term in (2- 

1 1). Note that the term has a varying near-wall behavior depending on the orientation of the 
rotating axis, but its lowest order is 0(y 2 ). Therefore, this term does not contribute to n* in 

Table 2.1, which only shows terms to 0{y) explicitly. 

From Table 2.1, one can see that £,y, Z),J and n* are the leading order terms in the near- 
wall region and £,y - n* is in balance with D up to O(y). To extend the high-Reynolds-number 
models to the near-wall region, appropriate expressions for n* and £ ( y are required so that they 

have the correct asymptotic behavior (Table 2. 1) in the near-wall region. Furthermore, according 
to Launder and Reynolds (1983), in the vicinity of the wall, the behavior of £,y / u^Uj has to satisfy 

the following kinematic relations 


2 2 

£j ] / U\ — £33 ! U 3 = £]3 / Wjl /3 = £ / k , 

(2-25a) 

£j 2 / mjm 2 = £ 23 / m 2 m 3 = 2 elk , 

(2-25b) 

£ 22 / u 2 = 4 £ / k 

(2-25c) 


These kinematic constraints are highly anisotropic, and they further indicate the difficulty of near- 
wall Reynolds-stress modeling, particularly for the 22 component. Another constraint for £,y is 

that its high-Reynolds-number model plus whatever near-wall corrections proposed should 
contract to 2e. 

A similar analysis carried out for (2-12) reveals that to O(y 0 ), Dj, D £ and y are the 
leading-order terms near the wall. Traditionally, D £ is neglected in the ad hoc modeling of the £- 

equation. The argument is that its effect can be accounted for in the model proposed for 
P l £ + Pj + P* - y in (2-12), i.e. the terms C £l f]Pe Ik- C e2 / 2 £ 2 Ik in (2-23). This reasoning is 

acceptable for high-Reynolds-number flows; but it is too restrictive for near-wall turbulence. The 
negligence of D £ is tantamount to the neglect of the pressure diffusion effect in the £ -equation, 
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which becomes very important in the near-wall region. Since the dissipation rate e influences the 
Reynolds stress u l u ] behavior in the computations through the interactions between £ tJ and UjUj , 

the importance of having a consistent near-wall model for the e-equation cannot be over 

emphasized. In view of this, near-wall modifications for the e-equation are necessary if proper 
near-wall modeling of is to be accomplished. 

Like other pressure-strain models invoking the high-Reynolds-number assumption, the SSG 
model is not asymptotically correct in the near-wall region. TheO(y°) term in the SSG model 
comes from the return term (Table 2.2) and is one order of magnitude larger than the exact term, 
which should be 0{y') (Table 1.1). Therefore, near-wall corrections are needed for the SSG 
model in order to make it applicable in the near-wall region. So are the dissipation rate tensor £ i} 

(2-22) and the dissipation rate equation developed for the high-Reynolds-number flows (2-23). In 
what follows, two near-wall Reynolds-stress models (NWRS and NNWRS) are presented in a 
parallel manner to highlight the difference between them. 

2.4. 1 NWRS Model With the above guideline, the NWRS model based on the SSG model 
was first developed by So et al. (1994). The pressure-strain correlation in the NWRS model can 
be written as 

n ij = -(C,e + C* P)b l} + C 2 e(b ik b kj -^n<5 y ) + C 5 k(b jm e ikm + b im e jkm )Q. k 

The last two terms are the near-wall corrections to the high-Reynolds-number SSG model (cf (2- 
26) to (2-21)). The Ilj- term, 

n Yj = (Qe + C[ P)bjj - C 2 e(b lk b kj - in<5, y ) + aiPij - ^P8 i} ) + 2 ykS i} , (2-27) 

where a* and /*are two new coefficients, is proposed to remove the 0(y°) terms in the SSG 
model (the Cj and C 2 terms in Table 2.2) and compensate for the insufficient anisotropy as a wall 
is approached. This term is multiplied by a damping function f w X to ensure that the near-wall 
correction diminishes away from the wall. The fl^ term, which is given by 
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where n,- = (0, 1, 0) is the wall unit vector, is included to improve the model predictions for the 

vv, uv and vw Reynolds-stress components in the near-wall region. In the course of validating 
the NWRS model, two different expressions for the damping function f w , are needed to predict 


correctly different types of flows with a wide range of Reynolds numbers. These expressions are 


fw, l = ex P 



(2-29a) 


fw, 1 = eX P 



(2-29b) 


where Re, = k 2 / (ve) is the local turbulent Reynolds number and 
A = \-{9l2)(b ij b ij -2b ij i^jk^ki) ' s anisotropy invariant. The choice of either (2-29a) or (2- 

29b) for the damping function depends on the flow Reynolds number and the type of flows 
considered. For example, when the model is applied to channel flow with Re T = 395, where 
Re T = u T h / v is the Reynolds number based on the friction velocity u r and the channel half width 
h, (2-29a) gives better results. On the other hand, when Re T is reduced to 180, (2-29b) gives 
better predictions. Overall, (2-29a) is suitable for most of the flows considered, except two cases 
of very low-Reynolds-number channel and pipe flows. The values of two additional near-wall 
constants a* and y* in (2-27) are specified to be a* = -0.29 and y* =0.065 to give the best 
predictions (compared with the DNS data of Kim et al. 1987 and Kim 1991) of the near-wall 
turbulence behavior for fully-developed channel flow at Re T = 180 and 395. 

A Similar approach is used to derive the near-wall correction function £ for the £ -equation 


(2-23). The improved £ -equation becomes 


n ^L--L( ill jL r - — :ii 

k dx k dx k . dx k , dx k s £ k ‘ dx i 


+ C.-P 


Q 2 f + C 

AC 


with 
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(2-30a) 


y , f ,£ s ,,£ 2 .,£fO 

f=/ "-V I T _a, tJ 1 

e = £ - 2vk / y 2 , (2-30b) 

£ = £-2v{d4k /<9y) 2 , (2-30c) 

where / lv2 = ex p[-fRe,/40j 2 is a damping function proposed to influence the decay of £ away 


from the wall, and L, M, and iV are model constants. Subsequent computations show that 
L = 2.25, M = 0.5 and A = 0.57 should be used. The modified dissipation rates £ and £ are 
introduced to ensure that the leading order terms in (2-30) have correct asymptotic behavior near a 
wall. 

The simple relation between the dissipation rate tensor e,y and the dissipation rate e for 


isotropic turbulence, (2-22), is no longer valid in the near-wall region. Further correction for (2- 
22) is required to account for turbulence anisotropy caused by the presence of the wall. Here, the 
relation given by Lai and So (1990), which satisfies the constraints (2-25) and contracts correctly 
to 2 e, is adopted: 


£ ij = f 1 ~ fw ,\ ) + fw,\ 


e UjUj +u i u k n k rij + UjU k n k n i + n l n j u k Ujn k n j 
k 1 + 3u k uin k ni / 2k 


(2-31) 


This relation consists of the isotropic part and the anisotropic correction part. The damping 
function f w ] is used to recover the Kolmogorov isotropic relation (2-22) at large Reynolds 

number. The near-wall asymptotic behavior of £,y given by (2-31) can be found in Table 2.3, 
which shows that (2-31) is only asymptotically correct to 0(y°) in the near-wall region. But each 
component has the correct leading order asymptotic behavior. 

2.4.2 NNWRS Model In the NWRS model presented in the previous section, all the near- 
wall corrections invoke wall dependent variables, such as found in (2-28), (2-30) and (2-31). As a 
result, the model is rather difficult to implement and may give incorrect predictions for flows with 
complex geometry. In this section, we develop a new near-wall Reynolds-stress (NNWRS) model 
with no wall normal dependence. 
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The pressure-strain term in the NNWRS model is modeled similarly as that in the NWRS 


model: 


n,y — (C|£ + C, P)by + C 2 £ ( b^ H ) + C^k(bj m e^ m + b uu c^ ni 

- «, - 1 ps,j ) - A < d,j - 1 ps,j ) - 2(r, + & n ' 1,2 )ks il + i ny. < 2 - 32a > 

where n* is again proposed to remove the incorrect lowest order O(y 0 ) terms in the SSG model 
and compensate for the insufficient anisotropy in the near-wall region, which is given by 

nJ=(C|£ + C,*W, y -C 2 e(b, k b kj , -in^ + aVi, - !*«(,•) + . (2-32b) 


and the damping function f w , is used to make Tly diminish in the high-Reynolds-number flow 


region. The damping function is proposed to have the following form, 

Av,i= ex p[— CV' 50 ) 2 - 


(2-33) 


by arguing that the near-wall modifications are due to the viscous effect and are needed when the 
local turbulence Reynolds number Re, is less than 150. This damping function is used for all 

types of flows. The new constants in (2-32) are chosen to be a = —0.36 and y = 0.072 to give 
good agreements with DNS data of fully developed channel flows at Re T = 180 and 395 (Kim et 
al. 1987; Kim 1991). Note that the wall-dependent near-wall correction term nj 2 (2-28) for the 
pressure- strain term in the NWRS model is not included in FI y in the NNWRS model (cf (2-32) 

with (2-26)). 

In the NWRS model, the modified dissipation rate e and e in the £ -equation depend on 
local normal wall distance y (2-30). To remove the wall dependence in the £-equation requires a 
complete re-examination of the derivation of the equation. Most recently, a new dissipation rate 
equation without wall dependence has been proposed by So et al. (1997) for two-equation model 
to account for wall effects. We re-examined the derivation of this equation thoroughly and decided 
to adopt it with an additional anisotropic coefficient in the turbulent diffusion term as the £- 
equation in the NNWRS model. The final form for the £ -equation in the NNWRS model is 
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u t 


de 


dx, dx 


V 


de 

dxi 


+ Cp-3 V 


kj 

£3 k v dx 


+ 4~( C e Z u ' u t *C c] jP- C c2 f e y 

dx k e dXj k k 


£ f d~Jk 


k J 


(2-34) 


where f £ = 1 - 0.22exp|-(Re,/ 6) 2 ] and f w 2 = exp[-(Re,/ 40) 2 ] are the damping functions, and 
C E =0.12, C £l =1.5, C e2 = 19 and C £3 = 2.95 are model coefficients. 

Finally, the relation for the dissipation rate tensor given by Shima (1995) is adopted in the 

NNWRS model, which is 
2 

k'ij — — ^ij £ 0 — f u,i ) 4" f w,\k'ij 4" £;/ 5 (2-3 5 a) 


where 


£ i} — —UjUj 


£,, = - 

y 2 


d , dujUj 


dx 


dxi 


UiUj d dk 

— ( v ) 

k dx k dx k 


(2-35b) 

(2-35c) 


The first term in (2-35a) is the isotropic part which recovers (2-22) when Reynolds number is 
large, the second term is the anisotropic correction in the near-wall region, and the last term is an 

additional near-wall correction which redistributes the dissipation rate among different Reynolds 
stress components. It can be shown that £„ given by (2-35) contracts to 2 £ and satisfies kinematic 

constraints (2-25) except the 22 component. The near-wall behavior of this dissipation rate tensor 
has been analyzed by applying (2-24) and is summarized in Table 2.3 in comparison with the 
corresponding components in the NWRS model (2-31). Table 2.3 shows that the asymptotic 
behavior of the dissipation rate tensor given by (2-35) matches the exact asymptotic behavior up to 
O(y'), whereas the one given by (2-31) in the NWRS model only up to <9(y°). The near-wall 
variations of different dissipation rate components given by (2-3 1 ) and (2-35) are compared with 
the DNS data (Kim et al. 1987), and the results are shown in Figs. 2.1a - 2. Id. In these figures, 
£," = VBjj / Uj are the non-dimensional dissipation rate components, and y + = yu r / v is the non- 
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dimensional wall distance, where u r is the friction velocity . These figures show that for the 
normal stresses, uu, vv, ww, the dissipation rate tensor relation in the NNWRS model gives 
better agreement with the DNS data than that in the NWRS model, but for the shear stress uv, 
both relations give about the same results. 

An alternative proposal has been put forward by Cho et al. (1995), where the dissipation rate 
tensor was assumed to be given by 


e ij ~ T $ij £ 0 - fwA ) + fw, l £ tj ' 


(2-36a) 


with 


£ V1 ’ rrr 
C IJ 


X + Wl~ f I U m U j - 

k k nmi i k Jm k 


1 | ^ £ pq U p U q 

2 £ k 


(2-36b) 


d4k d4k 


p — 2 v- 

c m« ~) -k 

dx t „ ox„ 


v m 


(2-36c) 


The asymptotic behavior of each component of this dissipation rate tensor and the corresponding 
overall behavior of the Reynolds-stress equations is given in Table 2.3 and 2.4 together with the 
NWRS and NNWRS models, respectively. According to these tables, this model (2-36) is only 
better than the NNWRS model in the 22 component. If (2-36) is used to replaced (2-35) in the 
NNWRS model, the resultant calculations yield better predictions in the near-wall region for the 
two channel flows considered. However, the results are not as good as those given by the NWRS 
model (Figs. 2.1 and 2.2). The improvement shows the importance of the 22 component in wall- 
bounded flows, thus the better performance of the NWRS model over the NNWRS model. 
Essentially, there is no difference in the overall predictions of the channel flows (Figs 2.2a - 2. 2d) 
when (2-35) and (2-36) are used. Therefore, the dissipation rate tensor (2-35) is adopted in the 
NNWRS model. 


2.5 Concluding Remarks 

In this section, we have presented two near-wall Reynolds-stress models, NWRS model and 
NNWRS model. The first one involves wall-dependent variables and the second one does not. 
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These near-wall models are developed from the high-Reynolds-number model by incorporating the 
near-wall modifications into the expression for the pressure-strain correlation, into the relation for 
the dissipation rate tensor in the Reynolds-stress equations, and into the modeled dissipation rate 
equation. Damping functions are used to ensure that the near-wall modifications will not affect the 

flow field in the region far away from the wall. The overall behavior of the Reynolds-stress 
equation in the near- wall region is determined by £,y - n* and is presented in Table 2.4 for both 

NNWRS and NWRS models. It shows that the modeled Reynolds-stress equations in both 
models match the exact equation to 0(y°). For those Reynolds-stress components (the 1 1, 33 
and 13) that are not correlated with the wall normal direction, they are better predicted by the 

NNWRS model than by the NWRS model. In contrast, the NWRS model, with additional near- 
wall correction term (see (2-28)), gives better prediction for those components (22, 12 and 

23) correlated with the wall normal direction. With no wall dependent variable in the model and 
only one expression for the damping function used in the pressure-strain formula, the NNWRS 
model is more general and flexible for flows with complex geometry than the NWRS model. The 
complete set of modeled equations and constants used in both models are given in Appendix B. 
The validations of the new developed near-wall Reynolds-stress model for different flow 
conditions will be presented in the following sections. 
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Table 2. 1 Near-wall asymptotic behavior of each term in Reynolds stress transport equation. 
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Efl 
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mm 

e v 
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n 

0{y 3 ) 

0(y 3 ) 

2vaf + \lva^a 2 y 
+ 0 (y 2 ) 

0(y 3 ) 

2va 2 + 8 va]a 2 y 

+ 0(y 2 ) 

-4 va ] a 2 y 

+0(y 2 ) 

22 

0(y 5 ) 

0(y 5 ) 

12 vb 2 y 2 

+0(y 3 ) 

0(y 5 ) 

8 vb]y 2 

+0(y 1 ) 

-4 vbly 2 
+ 0 (y 3 ) 

33 

0(y 2 ) 

0(y 3 ) 

2 VC 2 + 12 vc]C 2 y 

+ 0(y 2 ) 

0(y 3 ) 

2 vc 2 + 8 vc]C 2 >' 

+ 0(y 2 ) 

-4 vcjc 2 y 

+0(y 2 ) 

12 

0(y 4 ) 

0(y 4 ) 

6 va\b 2 y 

+0(y 2 ) 

0(y 4 ) 

4 Va,b 2 y 

+0(y 2 ) 

~2va^b 2 y 

+0(y 2 ) 

23 

0(y 4 ) 

0{y 4 ) 

6 vb 2 c^y 

+0(y 2 ) 

0(y 4 ) 

4 vb 2 c,y 

+0(y 2 ) 

-2 v/? 2 C] v 

+0(y 2 ) 

13 

0{ y 3 ) 

0(y 3 ) 

6 v(fl]C 2 + « 2 C 1 )y 

+ 0(y 2 ) 

0(y 3 ) 

2 va]C] 

+4 v(a^c 2 + a 2 C\ )y 

+ 0(y 2 ) 

- 2 v(a]C 2 + a 2 C] )y 

+ 0(y 2 ) 
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Table 2.2 Near-wali asymptotic behavior of each term in the SSG model for pressure-strain correlation FI ( y . 



24 






























Table 2.3 Near-wall asymptotic behavior of each component in the dissipation rate tensor e- . 
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NNWRS 
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Table 2.4 Near-wall asymptotic behavior of Ey n* (representing the overall behavior of the 
Reynolds-stress equation in the near-wall region). 
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NNWRS 
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Fig. 2. 1 Near-wall behavior of dissipation rate tensor £- compiled from DNS data for fully 
developed channel flow at Re T = 180 : (a) (b) £22 (c) £ 33 , and (d) £{ 2 . 
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3 Validation Against Simple Wall-Bounded Turbulent Flows 


3.1 Introduction 


In this section, the newly developed near-wall Reynolds stress (NNWRS) model is verified 
by comparing the model predictions with the DNS results and laboratory measurements for fully 
developed channel flow, pipe flow, and plane Couette flow, and for zero pressure gradient 
boundary layer flow over a flat plate. Further verification of the model predictions for complex 
turbulent flows are carried out in the following sections. For comparison, the corresponding 
predictions from the NWRS model are also presented in this section. The DNS experiments at 
low Reynolds numbers provide valuable data to validate turbulence models, especially in the near- 
wall region where measurement inaccuracy often arises. The simple flows considered in this 
section can be classified as internal and external flows. Channel flow, pipe flow and Couette flow 


are internal flows, whereas flat plate boundary-layer flow is an external flow. 

In a fully developed channel flow, the mean velocity is governed by 

1 dP d ( dU — ) 

= — v uv 

p dx dyy By j 


( 3 - 1 ) 


which shows that the constant pressure gradient balances the gradient of the total shear stress 
(including the viscous and turbulent stress), and acts as the driving force of the flow. The mean 
flow equation governing fully developed pipe flow is the same as (3-1) but written in cylindrical 
coordinate. The momentum equation for a fully developed Couette flow is governed by 
d f dU 



where the pressure gradient is zero, and the flow is driven by the moving wall. Consequently, the 
total shear stress in a Couette flow is constant in the entire flow domain. The pressure gradient in a 
boundary-layer flow over a flat plate is also zero, and equation (3-2) applies to the near-wall 
region. Away from the wall, the inertial force becomes important, and the momentum equation has 
the form 
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(3-3) 


u*L + v™ 

dx dy 


d( du 

V dy 


dy 


■ uv 


In wall-bounded flows, two regions can be identified; the inner layer and the outer layer (or 
the defect layer). The overlap of these two layers gives arise to the logarithmic layer where the log 
law of the wall applies (Fig. 3.1). Thus, the mean velocity in the log layer can be expressed as 
U + = K~ x In y + +B , (3-4) 

where U + = U / u T , y + = yu r / v are the non-dimensional mean velocity and wall normal 
distance, respectively; K is the von Karman constant, and B is another constant, with K~ 0.4 and 
5< 5 <5.5 depending on the type of flows. It is essential for a turbulence model to predict 
correctly the log-law region for both internal and external wall-bounded flows. Very close to the 
wall (y + < 5) is the viscous sublayer, where viscosity dominates and velocity varies linearly with 
the wall normal distance. In the outer layer (or defect layer), far away from the wall, the flow field 
is entirely turbulent, and fully developed channel, pipe and Couette flows show little wake 
characteristic compared with boundary layer flow with zero pressure gradient. 

In addition to the comparisons of the predicted mean and turbulent flow fields with existing 
data, the NNWRS model is also tested by examining its ability to reproduce Reynolds number 
effects on simple flows. Both the mean and the turbulent flow fields are influenced by the flow 
Reynolds number. So et al. (1996) demonstrated that the NWRS model is capable of assessing 
Reynolds number effects on internal and external flows. Following So et al. (1996), we shall also 
examine the ability of the NNWRS model to replicate Reynolds number effects in simple flows. 
The Reynolds number range considered varies from Re T = 180 to Re T = 8758 for fully-developed 
channel/pipe flows and from Re 0 = 1410 to Re e =2420 for flat plate boundary layers, where 
Re T = u T h / v, R v, h is the half-width of the channel or pipe radius, 6 is the 
momentum thickness and U ^ is the freestream mean velocity. 
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3.2 Fully-Developed Channel/Pipe and Couette Flows 

For fully developed channel/pipe and Couette flows, the modeled equations in NNWRS and 
NWRS models reduce to ordinary differential equations (see appendix C.1-C.3). Therefore, they 
are relatively easy to solve. Here an iterative scheme is used to solve the governing equations with 
the following convergence condition, 

<J)" +1 _ cp" 

J J 


where <X>” denotes a dependent variable at the jth grid point in the nth iteration. Only the lower half 

flow domain is considered in fully developed channel/pipe flows because the flow fields are 
symmetric about the centerline. The grid points are clustered close to the wall and stretched out 
away from the wall. For Couette flow, the entire flow domain is computed with specified mean 
velocity at the moving wall. The grid points are distributed symmetrically about the centerline with 
more points near the walls. The number of grid points used in channel/pine flow computations 
varies from 1 80 to 250 depending on the flow Reynolds number, and about twice the number of 
grid points are used in Couette flow computations. The first grid point away from the wall is 
located at y + = yu r / v = 1, where y is the wall distance. 

The no-slip boundary conditions are used to specify the mean velocity and turbulent stresses 
at a moving ( U m > 0) or stationary ( U m =0) wall, 

U = U m , uu = vv = ww — uv — 0 , (3-7a) 


< 10 " 


(3-6) 


and the following expression, 

£w= 2 v{djktdyf , 


(3-7b) 


is used to specify the dissipation rate £ at the wall. In fully developed channel/pipe flows, the 
boundary conditions at the centerline can be given as 
dU du 2 dv 2 


dy dy dy 


dw 2 de — _ 

^T“ MV = 0 


(3-7c) 


34 



due to the symmetry of the flow field (in case of pipe flow, cylindrical polar coordinates are used 
and (3-7c) with y + = R e T -r + , where r + is the radial coordinate, gives the corresponding 

boundary conditions at the centerline of the pipe ( r + = 0). 

3.2.1 Fully-developed channel flows For fully developed channel flow, two DNS cases 
studied by Kim et al. (1987) at Re T = 180 and Kim (1991) at Re r = 395 are used to test the 
models. Two damping functions, (2-29a) for Re T = 395 and (2-29b) for Re T = 180, are required 

for the NWRS model to give a reasonable agreement with the DNS results. On the other hand, 
only one damping function, (2-33), is used in the NNWRS model for both Re T = 180 and 395. 

The model predictions of the mean velocity profiles over the half channel width at Re r = 1 80 
and 395 are plotted in Figs. 3.2a - 3.2b. Both the NNWRS and NWRS model results agree very 
well with the DNS results (Kim et al. 1987; Kim 1991). Near the centerline, the velocity profiles 
predicted by the NNWRS model slightly deviate from those given by the DNS data. In the log 
layer, the von Karman constant (3-4) is determined with the procedure outlined by So et al. 
(1994). The NNWRS (NWRS) model gives K = 0.38 ( K = 0.40) at Re T = 180 and K = 0.40 
( K = 0.39) at Re T = 395 . Compared to K = 0.40 given by the DNS for both cases, the NNWRS 
slightly underpredicts the K value at Re T = 180, but gives same K value at Re T = 395. In 
contrast, the NWRS model underpredicts the K value at Re T = 395, but replicates the same DNS 
result at Re T = 180. Overall, the mean velocity profiles in fully developed channel flow given by 
both models, including the log-law behavior, are in good agreement with the DNS results. 

The comparisons between the model predictions for turbulence quantities in fully developed 
channel flow at Re r = 180 and Re T = 395 and the DNS results are presented in Figs. 3.3 - 3.6. 

Figs. 3.3a-3.3b show the turbulent kinetic energy profiles over the entire computational domain, 
i.e., half the channel width. The agreement between the NWRS model predictions and the DNS 
results is very good for both Reynolds numbers. The agreement between the NNWRS model 
results and the DNS results is reasonable, although the peak values of k + are underpredicted by 
the NNWRS model in both cases. Away from the wall, the NNWRS model gives slightly higher 
turbulent kinetic energy k + at Re T = 180, but slightly lower k + at Re r = 395. Figs. 3.4a-3.4b 
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plot the dissipation rate profiles for Re T = 180 and Re T =395 in the near-wall region, 
respectively. Both models essentially capture the variations of the dissipation rate predicted by 
DNS, e.g., a maximum value at the wall, a small plateau in a region close to the wall, and gradual 
attenuation away from the wall. The NNWRS and NWRS models overpredict the dissipation rate 
in the near-wall region, but agree well with the DNS away from the near-wall region. Figs. 3.5a- 
3.5b show the comparison of the predicted turbulent shear stress profiles with the DNS data. The 
shear stress uv + appears in the equation for mean flow (3-1), and therefore it is crucial for the 
correct prediction of the mean velocity. These figures show that the agreement between model 
predictions of the turbulent shear stress and the DNS results is excellent in both cases. Finally, 
Figs. 3.6a-3.6b plot the root-mean-square of the turbulent fluctuations (square root of the turbulent 
normal stresses) in the streamwise, wall normal and transverse directions, u+ ms , v^ ms , w+ ms , 

respectively. Compared to the DNS results, in the near-wall region, both models underpredict 

but overpredict v^ ms and w+ ms . The NWRS model predictions of the turbulent fluctuations 

are in better agreement with the DNS data than those from the NNWRS model. Consequently, the 
NWRS model predicts more accurately the anisotropy of the near-wall turbulence than the 
NNWRS model. Away from the wall, turbulence becomes more isotropic, and both models give 
almost the same results. 

According to asymptotic expansions (2-24), we have k + / e + = 0.5y +2 + 0(y +i ) in the near- 
wall region. To further examine the performance of the models in the near-wall region, we plot 
model predictions of k + / £ + versus y +2 in Fig. 3.7 to see if k + / (e + y +2 ) = 1 / 2 is true in the 
near-wall region. Fig. 3.7 shows that k + / (£ + y +2 ) = 1/2 is indeed true in the near-wall region for 
both models. 

3.2.2 Fully-developed pipe flows Fully developed pipe flow also has a constant pressure 
gradient in the streamwise direction. It is marginally more complicated than fully developed 
channel flow because of the cylindrical geometry. Four experimental data sets with a large range 
of Reynolds number are selected to test the models: they are the measurements of Durst et al. 
(1993) with Re r = 250, Schildknecht et al. (1979) with Re T = 489, and Laufer (1954) with 
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Re r = 1052 and 8758. In terms of bulk Reynolds number Re = U m D/ v, it ranges from 7,500 to 
500,000 in these four cases. Thus, the Reynolds number considered is quite high compared to 
many test cases investigated by other researchers. The number of grid points required for grid 
independent results are found be 108. In the NWRS model, damping function (2-29b) is used for 
the lowest Reynolds number case Re r = 250, and (2-29a) is used for the other three high 

Reynolds number cases. 

Figs. 3.8a - 3.8d show the comparison between the predicted and measured mean velocity 
profiles for the four cases considered. The NNWRS model gives an excellent prediction of U + in 
the near-wall region for all cases and in the entire flow region for Re T = 250 case, but it slightly 

underpredicts the velocity in the region close to the centerline as the Reynolds number increases. 
Even with the use of damping function (2-29b) for the low Reynolds number flow, the NWRS 
model still underpredicts the mean velocity in the log-law region for the Re r = 250 case. For the 

other three higher Reynolds number cases, the NWRS model also underpredicts the mean 
velocities in the region close the centerline, although its predictions are slightly better than those 
obtained from the NNWRS model. The predicted and measured von Karman constants for these 
four cases are presented in Table 3.1. Overall, both models reproduce well the mean velocity 
profiles including the log-law behavior in fully developed pipe flow. 

Figs. 3.9 - 3.12 compare the predicted turbulence quantities with available experimental data. 
In each figure, four panels are presented for the four different Reynolds number cases. In Figs. 
3.9a - 3.9d, the peak value of the turbulent kinetic energy k + is underestimated by the NNWRS 
model except for the Re r =250 case, whereas it is well predicted by the NWRS model in all 

cases. Away from the wall, both models give almost the same results, which agree with the 
measurements reasonably well. Both models give almost the same predictions of the dissipation 
rate (see Figs. 3.10a - 3.10d), except in the near-wall region where the dissipation rate predicted 
by the NWRS model has a more noticeable plateau. The predicted dissipation rate has a maximum 
value at the wall (Figs. 3. 10a - 3. lOd), which disagrees with the experimental data in the near-wall 
region (Figs. 3.10b - 3.10c). Away from the wall, the agreement between the model predictions 
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and measurements is very good. Figs. 3.1 la - 3.1 Id show the excellent agreement between the 
model predictions of the turbulent shear stress and the measurements whenever the measured data 
are available. As in channel flow, model predictions underestimate the turbulence fluctuation in 
the streamwise direction u+ ms in the near-wall region, whereas overestimate the turbulence 
fluctuations in the wall normal and transverse directions and w^ ms (Figs. 3.12a - 3.1 2d). The 
exception is the highest Reynolds number case Re T = 8758 in which the w+ ms is underpredicted 
by both models (Fig. 3.1 2d). Fig. 3.13 plots k + / £ + as a function of y +2 for the four cases with 
Re T ranging from 250 to 8758. As in channel flow, the predicted k + / e + in pipe flow follows 
the line representing k + / (e + y +2 ) = 1 / 2 in the near-wall region. 

3.2.3 Fully-developed plane Couette flows In plane Couette flows, one wall is stationary and 
the other is moving. A fully developed state is reached when the moving wall is dragged with a 
constant speed. The moving wall provides the energy to drive the flow. In the fully developed 
state, the total shear stress (viscous plus turbulent) is constant everywhere. Despite the seemingly 
simplicity of the flow, fully developed Couette flow is rather difficult to realize in laboratory, 
because of the difficulty in setting up the moving wall. But, it is a simple case for DNS (Lee and 
Kim 1991; Kristoffersen et al. 1993). The DNS data revealed some distinct features in Couette 
flow. For example, large scale eddies can be identified in the core region; flow in this region is 
quasi-homogeneous with turbulence production equal to dissipation rate but highly anisotropic. 
These features are quite different from those in Poiseuille flow driven by pressure gradient. 

Calculations are carried out for Couette flow with Reynolds number (based on half channel 
width) Re T = 170 and Re T = 625, and the results are compared with the DNS data given by Lee 
and Kim (1991) for the Re T = 170 case and the experimental data given by El Telbany and 
Reynolds (1980) for the Re T = 625 case. The damping function (2-29a) is found to be suitable 

for both cases, although the Reynolds number in Lee and Kim's (1991) case is rather low 
(Re T = 170). This indicates that plane Couette flow is less dependent on Re T than Poiseuille 

flow, which becomes more transparent by comparing the mean flow equations given in Appendix 
C. A term proportional to 1 / Re T appears in the non-dimensional mean flow equations for channel 
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(C-l) and pipe flows (C-14). However, no such term appears in the mean flow equation for 
Couette flow (C-27). 

Since all the turbulent quantities in Couette flow are symmetric about the centerline and 
(U / U m - 1 / 2) is anti-symmetric about the centerline (U m is the speed of the moving wall), we 

first present the model results in the lower half of the flow domain for both cases, and then show 
the results for the second case ( Re r = 625) in the entire domain. Figs. 3.14a-3.14b compare the 
model predictions of the mean velocity profiles with the DNS results for Re T = 170 and the 
experimental data for Re T = 625. Model predictions agree well with the data, although the 
NNWRS model slightly underpredicts the mean velocity in the core region. The predicted von 
Karman constant is very close to the K values from DNS and measurements (Table 3.2). 

In the near-wall region, the behavior of turbulent quantities is similar to their counterparts in 
Poiseuille flow (channel/pipe flows) (Figs. 3.15-3.18). Away from the near wall region, all the 
turbulent quantities rapidly becoming constant, and the turbulent flow field becomes homogeneous 
but remains isotropic (see Fig. 3.18). Again, the predicted k + / e + by both models in the near- 
wall region follow k + / (£ + y +2 ) = 1 / 2(Fig. 3.19). Fig. 3.20 compares the predicted U / U m , k + , 
e + , uv + , Up ms , v+ ms , w+ ms with the measurements for the Re T = 625 case over the entire flow 
domain. The mean velocity and the turbulent shear stress are accurately reproduced by both 
models. In the core region, both models underestimate the turbulent kinetic energy k + and the 
turbulent fluctuation in the wall normal direction . The turbulent fluctuation in the other two 

directions u+ ms and w+ ns are predicted quite accurately. 

3.3 Plane Boundary-Layer Flows with Zero Pressure Gradient 

For two dimensional boundary-layer flow, the governing equations of the models become 
parabolic (Appendix C.4). Wilcox's ( 1993) implicit, two-dimensional code developed for two- 
equation turbulence models is modified for the NNWRS and NWRS models. His code 
incorporates the compressible form of Levy-Lees transformation (Appendix C.5), which removes 
the singularity near the leading edge of the plate. In this transformation, dimensionless parameters 
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are introduced to represent the compressible effects. These parameters are set to zero and one for 
incompressible and isothermal flows, respectively. Thus, we still can use the same transformation 
to remove the singularity near the leading edge of the plate for incompressible turbulent boundary- 
layer flow. 

Computations are carried out by marching in the downstream (x) direction with iterations in 
the wall normal (y) direction. The step size in the downstream direction is adjusted according to 
how well the solution is converged at the previous cross-section location. In each cross-section, 
the distribution of the grid points follows a geometric progression ratio formula. After each 
iteration, the code checks for sufficient grid width. This is done by checking the difference of the 
mean velocity values between two successive grid points near the edge of the layer. If the 
difference is larger than about 10 -4 , a grid point is added to the domain. This ensures that all the 
wall normal derivatives will be zero near the boundary layer edge. Convergence is attained if the 
maximum errors in the mean and turbulence quantities between two successive iterations is less 
than 1CT 4 . 

The turbulence quantities at the edge of the boundary layer, uu e , vv e , ww e , uv e and e e , are 
determined by solving the following ordinary differential equations, 


U. 


' -C\£ e b x , + C 2 £ e {b} x + b\2 “ 3 _ J 


duu e 
€ ~lh 

dw e 

e !h 
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-C, 
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(3-8a) 

(3-8b) 

(3-Bc) 

(3-8d) 

(3-8e) 


which are reduced from the Reynolds stress equations in the models (Appendix C.4) under the 
assumption that all the wall normal derivatives vanish in the free stream. Therefore, 
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uu e , vv e , ww e , uv e , and e e given by (3-8) automatically satisfy the zero normal gradient 
conditions. 

At the initial marching location, a total of 150 grid points are used to cover the computational 
domain. As the boundary layer grows in the downstream, the number of grid points increases. In 
the numerical computations, the grids are distributed so that at least 15 points are located within 
y + <5, and more than 40 points are placed in the region 5<y + <65. Furthermore, the 
distribution ensures that the first grid away from the wall is placed within y + < 1 . 

Two boundary-layer flow cases are selected to test the models' ability to replicate simple 
external flows. Model results are compared with the DNS data given by Spalart (1988) at 
Re 0 = 1,410 and the detailed measurements of Karlsson and Johansson (1988) at Re e =2,420 
(Figs. 3.21 - 3.25). The predicted mean velocity profiles, which show the existence of a wake 
region in the outer layer, agree well with the DNS data and measurements (Figs. 3.21a - 3.2 1 b). 
The agreement between the predicted turbulent quantities k + , e + , uv + , u* ms , v^, and w+ ns and 
the available data for both cases is reasonable (Figs. 3.22 - 3.25). At Re e = 1,410, the turbulent 
quantities k + , uv + , u^ ms , v* ms , and w+ ms are first presented in the near-wall region and cross the 
boundary layer. In general, the NWRS model gives more satisfying predictions of the turbulence 
field in the near-wall region. The variations of the predicted turbulence quantities in the near-wall 
region are similar to those of the internal flows, with the exception that the maximum dissipation 
rate predicted by the NNWRS model is at a location very close to the wall, instead of the wall itself 
(Figs. 3.23a - 3.23b). As expected, all the turbulence quantities approach zero close to the edge of 
the boundary layer and match the freestream condition (Figs. 3.22 - 3.25). Fig. 3.26 shows that 
predicted ratio k + / e + in the near-wall region varies according to y +2 , as in the internal Hows. 

Note that unlike the internal flow computations in which Re T is specified, the boundary- 

layer computations are carried out by providing mean velocity and turbulent quantities at the initial 
location and free stream. Consequently, the wall shear stress, r w , is a predicted result rather than 
an input. The accuracy of this quantity or the skin friction coefficient Ce = 2x w / (pf/£) can be 

taken as an indication of the models' ability to predict external flows. Thus, in addition to 
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comparing the von Karman constant K, the calculated and measured C f are also compared (see 
Table 3.1). From Table 3.1, one can see that the predicted K and Cj agree well with the known 

results for both cases. 

3.4 Reynolds-Number Effects on Simple Turbulent Flows 

Using dimensional analysis, Millikan (1939) argued that Reynolds number effects could not 
be present in the mean velocity in wall-bounded turbulent flows. His argument is correct only in 
the near-wall region. Mellor and Gibson (1966) showed that Reynolds number has influence on 
the mean velocity in the outer region. Purtell et al (1981) pointed out that as Reynolds number 
decreases, the logarithmic region slowly disappears, while the viscous region remains unaffected 
by the decreasing Reynolds number. 

When local similarity arguments are applied to the higher-order turbulent statistics, they 
imply that, at least in the inner layer, individual second-order statistics at different Reynolds 
numbers would collapse into a single curve if they are nondimensionalized by using inner-layer 
scalings. Recently, direct numerical simulations, together with experimental measurements and 
analysis, have shown that Reynolds number effects on wall-bounded turbulent flows are evident, 
not only in the mean flow but also in the second-order statistics of the turbulence field. A most 
convincing demonstration of the Reynolds number effects on second-order statistics was given by 
Bandyopadhyay and Gad-el-Hak (1994), who showed that the location of the peak value of the 
turbulent shear stress, normalized by u/w T , increases with Reynolds number. So et al. (1996) 
have demonstrated that the NWRS model can predict the Reynolds number effects on mean and 
turbulence flow fields in simple flows with a wide range of Reynolds numbers. In what follows, 
we examine the ability of the NNWRS model to replicate the Reynolds number effects in wall- 
bounded simple turbulent flows. 

3.4.1 Reynolds -number effects on mean flow To show the Reynolds number effects on the 
mean flow, the predicted normalized mean velocity by the NNWRS model for internal (channel 
and pipe) and external (boundary layer) flows with different Reynolds numbers are plotted together 
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in Figs. 3.27a and 3.27b, respectively. Note that Couette flow is excluded in this section because 
it is less dependent on Reynolds number compared with the other internal flows. For comparison, 
the available DNS and experimental data are also plotted in Fig. 3.27. For channel/pipe flow, 
model predictions and the available data show little wake component, and the log layer extends all 
the way to the centerline (Fig. 3.27a). The extent of the log layer increases as Reynolds number 
increases. In the viscous layer, the data and model results are all very well correlated by a single 
curve U + -y + . On the other hand, for boundary layer flow, the model results and the available 
data show a wake component in the defect layer and a reduction of the log-law region as Reynolds 
number decreases (Fig. 3.27b). However, the mean velocity in the viscous layer is unaffected by 
Reynolds number and again is well approximated by U + = y + . 

In the log layer, the velocity profile varies according to (3-4). The von Karman constant K 
obtained by following the procedure outlined in So et al. (1994) for each individual case has been 
given in Table 3.1. From Table 3.1, one can see that Reynolds number has no influence on the 
von Karman constant in wall-bounded flows. In fact, the velocity profiles for different Reynolds 
numbers in the log layer can be well correlated by (3-4) with K = 0.40 and 5 = 5.2 for 
channel/pipe flow (Fig. 3.27a), except for the Re T =8758 case where the mean velocity is 
underpredicted by the NNWRS model, and with K = 0.41 and B = 5.2 for boundary layer flow 
(Fig 3.27b). Both models do a fair job of reproducing the universal inner-layer behavior and the 
prediction of K is within the error margin of its determination from experiments. In general, a 
value of 0.40 is obtained and this is in agreement with experimental and DNS data. 

For boundary-layer flow, the effects of Reynolds number on the shape factor H = 8* / 0 
(where 8* and 8 are the displacement and momentum thickness of the boundary layer, 

'y 

respectively) and skin friction coefficient Cj = 2r w / pU ^ are also examined. The predicted and 
measured H and Cj are reported in Table 3.1. The DNS and experimental data show that H and 
Cj decrease as Reynolds number increases, which is also predicted by the models. Both models 
predict the shape factor well. But the NWRS model gives a more accurate prediction of the skin 
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friction coefficient than the NNWRS model. Overall, Reynolds number effects on the mean flow 
are fairly well reproduced by the NNWRS model, as well as by the NWRS model. 

3.4.2 Reynolds -number effects on turbulence statistics Momentum transport in two- 
dimensional wall-bounded turbulent flows is mainly carried out by the Reynolds shear stress. 
Therefore, it is of paramount importance to model the shear stress behavior correctly. In fully 
developed channel/pipe flows, the normalized shear stress can be expressed as 

-w + =(\-y + /Re r )-dU + /dy + . (3-9) 

Assuming a universal velocity profile, it can be seen from (3-9) that the Reynolds number 
dependence of the shear stress in the inner layer is rather strong at low Reynolds number. Wei and 
Willmarth (1989) examined channel/pipe flows and found that the normalized shear stress at 
different Reynolds numbers does not collapse in the outer layer and the separation of the different 
profiles is still discernible in the inner layer. The calculated shear stress profiles for different 
Reynolds numbers are plotted together in Figs. 3.28a and 3.28b for channel/pipe and boundary- 
layer flows, respectively. Whenever the DNS and experimental data are available, they are also 
plotted in the figures for comparison. Fig 3.28 shows that in the logarithmic and outer regions, the 
shear stress profiles spread out as Reynolds number increases. In the viscous layer, the shear 
stress profiles at different Reynolds number in channel/pipe flow do not quite collapse into a single 
curve, whereas they do collapse into a single curve in boundary -layer flow. Thus, the influence of 
Reynolds number on the shear stress is more pronounced in channel/pipe flow than in boundary- 
layer flow. 

Another test of the ability of the NNWRS model to reproduce Reynolds number effects is to 
plot the location of the peak shear stress versus Re r ; for boundary layer flow, the reduced 
Reynolds number defined as Re* = u r S / v is used instead (Bandyopadhyay and Gad-el-Hak 
1994). The location of the peak shear stress, denoted by y p , can be determined from the shear 
stress profiles shown in Fig. 3.28. A log-log plot of y p versus Re r is shown in Fig. 3.29 for all 
the flows examined. According to Sreenivasan (1988), a linear relation should exist between In y p 
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and lnRe T . Bandyopadhyay and Gad-el-Hak (1994) showed that the straight line y p = 2(Re r ) 1/2 
correlates well with the experimental data they examined in a log-log plot. Fig. 3.29 shows that 
the calculated results by the NNWRS model also correlate well with the same straight line. 

The normalized production of kinetic energy, P + = -uv + (dU + / dy + ), is shown to be 
relatively independent of Reynolds number in the outer region of channel flow but not so in the 
inner region (Wei and Willmarth 1989). Furthermore, the location of the peak value of P + is also 
relatively independent of Reynolds number (Bandyopadhyay and Gad-el-Hak 1994), although the 
location of the peak shear stress varies with Re r (see Fig. 3.29). These conclusions can also be 
deduced from (3-9) and the definition of P + . As Re T — » P + reaches its maximum 1/4. The 
calculations and measurements of P + plotted versus In y + are shown in Figs. 3.30a - 3.30b. As 
expected, there are no Reynolds number effects on the calculated P + in the outer layer for 
channel/pipe flow and boundary-layer flow. In the inner layer, the dependence of P + on 
Reynolds number in internal and external flows is different (cf. Fig. 3.30a with 3.30b). In 
channel/pipe flow, P + , particularly its peak value, distinctly depends on Re T . Although the 
calculated location of the peak shear stress uv + varies with Re T (Fig. 3.29), there is little variation 
in the location of the peak value of the predicted P + . In both channel/pipe flow and boundary- 
layer flow, P + peaks at about y + = 10, which is consistent with y + = 12 given by 
Bandyopadhyay and Gad-el-Hak (1994) who examined different sets of data with a different 
Reynolds number range. The peak value of the calculated P + is close to 1/4 and approaching 1/4 
only at high Reynolds number. The Reynolds number also has more effects on the peak value of 
the shear stress than the peak value of the production of kinetic energy (cf. Fig. 3.30 with 3.28). 
The ability of the NNWRS model to predict these important features of wall-bounded flows is 
another indication of the validity of the model. 

The predicted turbulent dissipation rate £ + and the viscous diffusion of turbulent kinetic 
energy, D ^ = d 2 k + / dy +2 , are compared with DNS data and measurements in Figs. 3.31a - 
3.31b. Very near the wall, the dissipation rate is approximately balanced by the viscous diffusion. 
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According to the asymptotic analysis in Section 2, this suggests that the pressure diffusion is 
much weaker than the viscous diffusion and the dissipation rate in the viscous layer. The predicted 
viscous diffusion is not sensitive to the variation of Reynolds number, neither is the dissipation 
rate except at very low Reynolds numbers. 

3.5 Concluding Remarks 

In this section, the NNWRS model has been verified by comparing model predictions with 
the DNS data and measurements for fully-developed wall-bounded simple flows: channel flow, 
pipe flow, and Couette flow, and for boundary-layer flow. The results show that the NNWRS 
model predicts the mean and turbulent flow fields of internal and external flows reasonably well. 
The model is also able to capture the log-law region with a correct von Karman constant and 
replicate the Reynolds number effects on the mean flow and second-order statistics. The main 
deficiency in the model is that it predicts less degree of turbulence anisotropy in the near-wall 
region compared to the NWRS model. This may be attributed to the more general formulation in 
the NNWRS model, i.e., no wall normal has been used in the model. In the next section, we shall 
apply the NNWRS to study wall-bounded complex turbulent flows. 
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Table 3.2 Comparison of the predicted K with data in Couette flow. 


Data Source 

Re T 

Data 

NNWRS 

NWRS 

Lee and Kim (1991) 
(DNS) 

170 

0.40 

0.40 

0.39 

El Telbany and Reynolds 
(1980) (EXP) 

625 

0.39 

0.40 

0.40 
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Fig. 3.1 Sketch of different flow regions in a wall-bounded turbulent flow. 
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Fig. 3.2a Comparison between the predicted and DNS mean velocity profiles in the fully 

developed channel flow at Re T = 180. 



Fig. 3.2b Comparison between the predicted and DNS mean velocity profiles in the fully 

developed channel flow at Re T = 395. 
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Fig. 3.3a Comparison between the predicted and DNS turbulent kinetic energy in the fully 

developed channel flow at Re r = 180. 



Fig. 3.3b Comparison between the predicted and DNS turbulent kinetic energy in the fully 

developed channel flow at Re r = 395. 
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Fig. 3.4a Comparison between the predicted and DNS dissipation rate in the near-wall region of 

fully developed channel flow at Re T = 180. 



Fig. 3.4b Comparison between the predicted and DNS dissipation rate in the near-wall region of 

fully developed channel flow at Re T = 395. 
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Fig. 3.5a Comparison between the predicted and DNS turbulent shear stress in the fully developed 

channel flow at Re T = 180. 



Fig. 3.5b Comparison between the predicted and DNS turbulent shear stress in the fully developed 

channel flow at Re T = 395 . 
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Fig. 3.6a Comparison between the predicted and DNS turbulent fluctuations in the fully developed 

channel flow at Re T = 180. 



Fig. 3.6b Comparison between the predicted and DNS turbulent fluctuations in the fully developed 

channel flow at Re f = 395 . 
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Fig. 3.7 Predicted k + / e + vs y +2 in the near-wall region for fully developed channel flows. 
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Fig. 3.8a Comparison between the predicted and measured mean velocity profiles in the fully 

developed pipe flow at Re T = 250. 



Fig. 3.8b Comparison between the predicted and measured mean velocity profiles in the fully 

developed pipe flow at Re T = 489. 
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Fig. 3.8c Comparison between the predicted and measured mean velocity profiles in the fully 

developed pipe flow at Re r = 1052. 



Fig. 3.8d Comparison between the predicted and measured mean velocity profiles in the fully 

developed pipe flow at Re T = 8758. 
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Fig. 3.9a Comparison between the predicted and measured turbulent kinetic energy in the fully 

developed pipe flow at Re T = 250. 



Fig. 3.9b Comparison between the predicted and measured turbulent kinetic energy in the fully 

developed pipe flow at Re r = 489. 
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Fig. 3.9c Comparison between the predicted and measured turbulent kinetic energy in the fully 

developed pipe flow at Re T = 1052. 



Fig. 3.9d Comparison between the predicted and measured turbulent kinetic energy in the fully 

developed pipe flow at Re T = 8758. 
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Fig. 3. 10a Comparison of the predicted dissipation rate in the near-wall region of fully developed 

pipe flow at Re r = 250. 



Fig. 3.10b Comparison between the predicted and measured dissipation rate in the near-wall region 

of fully developed pipe flow at Re T = 489. 
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Fig. 3.10c Comparison between the predicted and measured dissipation rate in the near-wall region 

of fully developed pipe flow at Re r = 1052. 



Fig. 3. lOd Comparison of the predicted dissipation rate in the near-wall region of fully developed 

pipe flow at Re T = 8758. 
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Fig. 3. 1 la Comparison of the predicted turbulent shear stress in the fully developed pipe flow at 

Re r = 250. 



Fig. 3.1 lb Comparison between the predicted and measured turbulent shear stress in the fully 

developed pipe flow at Re T = 489. 
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Fig. 3.1 lc Comparison between the predicted and measured turbulent shear stress in the fully 

developed pipe flow at Re r = 1052. 



Fig. 3.1 Id Comparison between the predicted and measured turbulent shear stress in the fully 

developed pipe flow at Re r = 8758. 
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Fig. 3.12a Comparison between the predicted and measured turbulent fluctuations in the fully 

developed pipe flow at Re r = 250. 



Fig. 3.12b Comparison between the predicted and measured turbulent fluctuations in the fully 

developed pipe flow at Re T = 489. 
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Fig. 3.12c Comparison between the predicted and measured turbulent fluctuations in the fully 

developed pipe flow at Re r = 1052. 



Fig. 3.12d Comparison between the predicted and measured turbulent fluctuations in the fully 

developed pipe flow at Re T = 8758. 


65 




y + 2 =(Re t -r + ) 2 


Fig. 3.13 Predicted k + / e + vs. y +2 = (Re T - r + ) 2 in the near-wall region for fully developed pipe 

flow. 
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Fig. 3.14a Comparison between the predicted and DNS mean velocity profiles in the fully 

developed Couette flow at Re T = 170. 



Fig. 3.14b Comparison between the predicted and measured mean velocity profiles in the fully 

developed Couette flow at Re T = 625. 
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Fig. 3.15a Comparison between the predicted and DNS turbulent kinetic energy in the fully 

developed Couette flow at Re T = 170. 



Fig. 3.15b Comparison between the predicted and measured turbulent kinetic energy in the fully 

developed Couette flow at Re T = 625 . 
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Fig. 3.16a Comparison of the predicted dissipation rate in the near-wall region of fully developed 

Couette flow at Re T = 170. 



Fig. 3.16b Comparison of the predicted dissipation rate in the near-wall region of fully developed 

Couette flow at Re r = 625. 
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Fig. 3.17a Comparison between the predicted and DNS turbulent shear stress in the fully 

developed Couette flow at Re T = 170. 



y + 


Fig. 3.17b Comparison between the predicted and measured turbulent shear stress in the fully 

developed Couette flow at Re T = 625. 
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Fig. 3.18a Comparison between the predicted and DNS turbulent fluctuations in the fully 

developed Couette flow at Re T = 170. 



Fig. 3.18b Comparison between the predicted and measured turbulent fluctuations in the fully 

developed Couette flow at Re r = 625. 
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Fig. 19 Predicted k + / e + vs. y +1 in the near-wall region for fully developed Couette flows. 
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Fig. 3.20a Comparison between the predicted and measured mean velocity profiles over the entire 

fully developed Couette flow region at Re T = 625. 


5 
4 
3 

k + 

2 

1 

0 

0 0.4 0.8 „ 1.2 1.6 2 

y/h 


Fig. 3.20b Comparison between the predicted and measured turbulent kinetic energy over the 
entire fully developed Couette flow region at Re r = 625 . 
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Fig. 3.20c Comparison of the predicted dissipation rate over the entire fully developed Couette 

flow region at Re T = 625. 
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Fig.3.20d Comparison between the predicted and measured shear stress over the entire fully 

developed Couette flow region at Re T = 625 . 
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Fig. 3.20e Comparison between the predicted and measured turbulent fluctuations over the entire 
fully developed Couette flow region at Re T = 625. 
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Fig. 3.21a Comparison between the predicted and DNS mean velocity profiles of boundary-layer 

flow at Re e = 1410. 



Fig. 3.21b Comparison between the predicted and measured mean velocity profiles of boundary- 

layer flow at Re 0 = 2420. 
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Fig. 3.22a Comparison between the predicted and DNS turbulent kinetic energy in the near-wall 

region of boundary-layer flow at Re 0 = 1410. 



Fig. 3.22b Comparison between the predicted and DNS turbulent kinetic energy cross the 

boundary layer at Re e = 1410. 
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Fig. 3.23 Comparison between the predicted and measured turbulent kinetic energy cross the 

boundary layer at Re e = 2420. 
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Fig. 3.24a Comparison between the predicted and DNS dissipation rate in the near-wall region of 

boundary-layer flow at Re 0 = 1410. 



Fig. 3.24b Comparison between the predicted and measured dissipation rate in the near-wall region 

of boundary-layer flow at Re e = 2420. 
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Fig. 3.25a Comparison between the predicted and DNS turbulent shear stress in in the near-wall 

region of boundary-layer flow at Re 0 = 1410. 



Fig. 3.25b Comparison between the predicted and DNS turbulent shear stress cross boundary 

layer at Re 0 = 1410. 
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Fig. 3.26 Comparison between the predicted and measured turbulent shear stress cross boundary 

layer at Re 0 = 2420 . 
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Fig. 3.27a Comparison between the predicted and DNS turbulent fluctuations in the near-wall 

region of boundary layer flow at Re 0 = 1410. 



Fig. 3.27b Comparison between the predicted and DNS turbulent fluctuations cross boundary 

layer at Re e = 1410. 
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Fig. 3.28 Comparison between the predicted and measured turbulent fluctuations cross boundary 

layer at Re 0 = 2420. 



Fig. 3.29 Predicted k + / e + vs. y +2 in the near-wall region for boundary layer flows. 
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Fig. 3.30a Comparison of calculated mean velocity plotted in inner-layer variables to show 
Reynolds number effects in channel and pipe flows. 



Fig. 3.30b Comparison of calculated mean velocity plotted in inner-layer variables to show 
Reynolds number effects in boundary layer flows. 
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Fig. 3.31a Comparison of calculated turbulent shear stress «v + with data to show Reynolds 

number effects in channel and pipe flows. 
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Fig. 3.31b Comparison of calculated turbulent shear stress uv + with data to show Reynolds 

number effects in boundary-layer flows. 
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Fig. 3.33a Comparison of the calculated production of turbulent kinetic energy P + with data to 
show Reynolds number effects in channel and pipe flows. 
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Fig. 3.33b Comparison of the calculated production of turbulent kinetic energy P + with data to 
show Reynolds number effects in boundary layer flows. 


89 



Calculation 

Data 

Re 

X 


o 

180 



250 


A 

395 




489 

- — - 

□ 

1052 

8758 



Fig. 3.34a Comparison of the calculated turbulent dissipation rate e + and diffusion D 
turbulent kinetic energy with data to show Reynolds number effects in channel and pipe flows. 
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Fig. 3.34b Comparison of the calculated turbulent dissipation rate £ + and diffusion D 
turbulent kinetic energy with data to show Reynolds number effects in boundary-layer flows. 
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4 VALIDATION AGAINST TURBULENT FLOWS WITH SWIRL AND 

SYSTEM ROTATION 


4.1 Introduction 

Swirling and rotating turbulent flows are often encountered in engineering applications. In 
swirling flows, the induced swirl velocity component changes the paths of fluid particles from 
parallel to spiral. The flow pattern is affected by the centrifugal force associated with local 
streamline curvature. Recirculation may appear in the central region of the pipe close to the inlet. 
In rotating flows, the rotation effect on flow patterns depends on the orientation of the rotating 
axis to the mean flow plane, i.e. the relative importance of the centrifugal force to the Coriolis 
force. Flow in a rotating pipe is mainly subject to the centrifugal force because the rotating axis is 
parallel to the mean flow direction. On the other hand, only Coriolis force is dynamically 
important to fully developed channel flow with rotating axis perpendicular to the mean flow plane. 
Therefore, with the presence of swirl and rotation, the flows become much more complicated than 
those discussed in Section 3. Although the present wall-independent near-wall Reynolds stress 
(NNWRS) model is developed for flows with complex geometry, the model is also expected to 
perform well for swirling and rotating flows. In this section, its ability to replicate such complex 
flows is tested and compared with that of the NWRS model. 

4.2 Turbulent Swirling Flows in a Straight Pipe 

4.2. 1 Background Earlier studies on swirling flows were mainly concerned with the mean flow 
and pressure drop measurements and the decay of swirl along the pipe (Kreith and Sonju 1965; 
Backshall and Landis 1969; Yajnik and Subbaiah 1973; Murakami et al. 1976; Padmanabhan and 
Janek 1980; Ito et al. 1980; Kito 1984; Kito and Kato 1984). Little attention was paid to the decay 
of the turbulence field. Measurements on the evolution of the turbulence field were carried out 
only recently (Algifri et al. 1987; Kitoh 1991; Parchen et al. 1993). 

Swirling flows in a straight pipe can be generated by rotating an inlet section of the pipe 
(Weske and Sturov 1974) or by a vane swirler installed at the entrance of the pipe (Kitoh 1991). 
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The swirl intensity can be measured by the swirl number S defined as the ratio of the angular 
momentum to the axial momentum. For axisymmetric flows, swirl number can be expressed as 
(Kitoh, 1991) 


2 f° UWr 2 dr 

c _ 

O — in 


r 2 I ! 2 

r 0 U m 


(4-1) 


where U and W are the mean velocity components in the axial and azimuthal, respectively; U m is 
the bulk mean axial velocity; r and r 0 are the radial position and the pipe radius. 

Fig. 4.1 shows a typical azimuthal velocity profile W in axisymmetric swirling pipe flows. 
In the so called forced-vortex region, the velocity profile resembles that associated with rigid-body 
rotation, i.e., the velocity linearly increases with increasing distance from the centerline. The mean 
vorticity in the axial direction is a constant in this region. Outside the forced-vortex region, the 
velocity starts decreasing and matches the wall boundary condition through the boundary layer. 
The region between the boundary layer and the forced-vortex region is called the free-vortex region 
because the axial vorticity strength is almost zero there. Table 4.1 summaries the characteristics of 
swirling flows in forced- and free-vortex region. In forced vortex region, the rotational strain 
associated with swirl (dW / dr - W / r) / 2 is zero and swirl does not provide additional turbulent 
production. On the other hand, in the free-vortex region, the rotational strain is not zero and swirl 
provides extra production. Turbulence is thus stabilized in the forced- vortex region and 
destabilized in the free vortex generated by swirl. The extent of the forced- and free-vortex 
region in swirling flows depends on swirl intensity and the way swirl is generated. For instance, 
the swirling flows in Weske and Sturov's (1974) experiments , generated by a rotating gird in the 
inlet section, are forced-vortex-dominated, whereas the swirling flows generated by a vane swirler 
in Kitoh's ( 1 99 1 ) experiments are free-vortex-type. 

Flow characteristics near the pipe centerline also depends on the swirl number S. For flows 
with S < 1 , a reverse flow region usually is not observed because the pressure depression is not 
strong enough to create a reverse flow. For flows with S > 1 , a reverse flow region may appear 
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depending on whether other conditions, such as Reynolds number, inlet conditions and the way 
swirl is generated, are favorable to the formation of such region or not. For example, the 
experiments of Kitoh (1991) showed a reverse flow region existed, whereas the experiments of 
Weske and Sturov (1974) did not indicate such region existed, although the initial swirl number S 0 
is greater than one in both cases. 

Swirl decays along the pipe as a result of wall friction, and the mean azimuthal velocity 
profile alters as flow approaches to the fully-developed state downstream. Swirl decay rate was 
found to be dependent on inlet swirl intensity, Reynolds number, and pipe roughness (Seno and 
Nagata 1972; Baker 1967; Padmanabhan and Janek 1980). For small swirl number, Kitoh (1991) 
derived an exponential formula to predict the attenuation of axisymmetric swirling flow along the 
pipe. For large swirl number, swirl decay rate can be obtained only by experiments or numerical 
computations. 

The unique features of swirling flows have already been made use of in a number of 
engineering applications. Prominent among them are flames in gas turbines and furnace 
combustors in which swirl strongly contributes to efficiency of combustion by enhancing mixing, 
and to flame stability through recirculation. The primary mixing enhancement is attributed to 
higher levels of turbulence generated by the additional mean shear strain. Table 4.2 compares the 
mean shear strain components in swirling flows and parallel flows. 

Turbulence in swirling flows is anisotropic. The anisotropy results from the uneven 
weighting of swirl effect on turbulent normal and shear stresses. The degree of anisotropy 
depends on the swirl intensity and the location in the flow. A rather careful analysis of this 
anisotropic behavior was carried out by Lilley and Chigier (1971) using the mean flow 
measurements of the swirling free jet experiment of Lilley and Chigier (1967). They found that 
depending on swirl number, the ratio x n /x rB varies from 2 to more than 8, where t rx and x r6 
are the rx and rd components of the turbulent shear stress, respectively. Consequently, only 
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those more advanced turbulence models that can replicate stress anisotropic behavior are expected 
to give satisfactory predictions of swirling flows. 

So far, most modeling work on swirling flows have been carried out by using modified 
versions of two-equation models (Sloan et al. 1986). With the isotropic turbulence assumption, 
two-equation models need to be modified before they are used to predict anisotropic swirling 
flows. Great efforts have been put on the improvement of swirl effect in two-equation modeling. 
The improvement was often marginal and obtained in a single tested case without further validation 
in others. This necessitates the need of using other turbulence modeling beyond the Boussinesq 
approximation. In other words, Reynolds stress modeling is needed for a more accurate prediction 
of swirling flows. A small number of attempts have been made to predict strongly swirling flows 
(free or confined) by using Reynolds-stress models (Jones and Pascau 1989; Hogg and Leschziner 
1989; Fu et al. 1988). The results show clearly the superiority of the Reynolds stress models to 
the modified two-equation models. 

4.2.2 Boundary conditions and numerical implementation For axisymmetric swirling flows 
in a straight pipe, the model governing equations are elliptic. To solve the equations in the (x, r)- 
plane, boundary conditions at the inlet and outlet of the computational domain, along the centerline 
of the pipe and at the wall are required. The inlet boundary conditions for numerical computations 
are provided by the available experimental data at the first measured location (with interpolations if 
necessary). Turbulence kinetic energy k in at the inlet of the computational domain can be specified 

from the available turbulent normal stresses by 

, UUin ~t~ Win ~h WWin , A 

t = . (4-z) 

tn ^ 


The dissipation rate at the inlet is estimated from the turbulent kinetic energy according to 

i.3/2 

A 




xi 


(4-3) 


where / is the characteristic length scale of the case considered (such as the radius of the pipe in a 
swirling pipe flow), and X =0.02 is a parameter. The turbulent shear stresses at the inlet are 
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specified from the mean shear strain rates through the Boussinesq approximation. At the outlet of 
the computational domain, negligible diffusion condition in the axial direction is found to be the 
most appropriate outlet boundary condition. This implies that the second derivatives of the 
dependent variables can be assumed to be zero for a sufficient long computational domain (100 
diameters in the present study). 


Along the centerline of the pipe, zero normal gradient is specified for all dependent variables 


except those that they are zero themselves, 

d(U, W, mm, vv, ww, mw, e) 
=0 , an 

dr 


= 0 , and V = mv = vw = 0 



The TEACH code for incompressible flows by Gosman and Ideriah (1976) is adopted for the 
numerical computations. This code uses the finite volume (cell) method to discretize the transport 
equations and adopts the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) 
algorithm of Patankar and Spalding (1972) to solve the discretized system of equations. A detailed 
account of the finite volume method and SIMPLE algorithm can be found in Patankar (1980), 
among many other references. Here, a brief description is given to the finite volume method and 
the solution procedure. 

The time-averaged transport equations of mean and turbulence quantities can be rewritten in 
the following form 

V»(U 0 )=V*(r <p V 0 ) + S 0 , (4-6) 

where d> is a dependent variable representing £/ ( , e, and m,m 7 ; is the effective diffusion 
coefficient, and S<j, represents the rest of terms that cannot be included into convection and 
diffusion terms. The flow domain is divided into a series of control volumes according to the 
given grid points. Equation (4-6) is integrated within each control volume and the result is 
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expressed in terms of the unknown variable at neighboring grid points. In a control volume 
centered at grid point P (Fig. 4.2), which is surrounded by the neighboring grid points E, W, N 
and S, the resulting discretized equation for O can be written as 

Ap<P p = A e O e + A W <P W + A n <P n + A S 0 S + SAV , (4-7) 

where the subscript of O indicates the evaluated location; coefficients As involve U i and r 0 , and 

depend on the method used to discretize the integrals corresponding to the convection and diffusion 
terms in (4-6); S is the average of the source term S t j> in (4-6) over the control volume, and AV is 
the volume of the cell. The averaged source term S is expressed formally as a linear function of 

O, 

S = S c +Sp^ p , (4-8) 

where S p is the coefficient of O,,, and S c is the part in 5 that does not depend on <J> P (formally). 
For the turbulent transport equations governing the Reynolds stresses and dissipation rate e , the 
averaged source term S is a strong function of the dependent variable <t». One can come up with 
different expressions for S p and S c , especially for the complex source terms in the Reynolds- 
stress equations. It suffices to say that the way S is partitioned into S c and S P <t> P plays a very 
crucial role in finding the solution. In general, it is desirable to have a negative S P since a positive 
S P could cause divergence. 

In the present modeling of axisymmetric swirling flows, the hybrid scheme (a combination of 
the central difference and upwind schemes) described by Patankar (1980) is used. The control 
volume shown in Fig. 4.2 is used for all the dependent variables except the mean velocity in the 
axial and radial directions U and V. The control volumes for U and V are staggered as shown in 
Fig. 4.3a and Fig. 4.3b, respectively. 

Since the system of equations (4-6) are nonlinear with coefficients and source terms as 
functions of <5 itself and other dependent variables, the final solution is obtained by iteration. 
Under-relaxation iteration technique (Patankar 1980) is used to solve the system of equations. 
With the introduction of the under-relaxation parameter a, (4-7) can be rewritten as, 
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(4-11) 


— <P p =A e 0 e + A w 4> w + A„4>„ + A S 0 S + SAV + f L . a > A ^ 
a a 

where Op stands for the value of Op from the previous iteration. 

A tridiagonal matrix algorithm is used to solve the system of equations (4-1 1). A line-by- 
line iterative approach is adopted in alternating directions within the consecutive stages of the 
SIMPLE algorithm as described by Patankar (1980). The details of handling the mean velocity and 
pressure variables and associated difficulties can also be found in Patankar (1980). 

The implementation of the above numerical procedure to solve a Reynolds stress model is not 
a trivial task. For axisymmetric swirling flows, we still have eleven equations for eleven 
unknowns: three mean velocity components, six Reynolds stresses, one dissipation rate and one 
mean pressure variable. Obtaining the solutions of these highly nonlinear and coupled transport 
equations requires a great number of iterations. Intermediate results of a near-wall two-equation 
model are used to initialize the iterations, which usually lead to the converged solution faster. 
Converged solutions are obtained when the maximum residuals of the mass and momentum 
equations in the entire computational domain are less than a small number 8 = 0( 10~ 3 ), i.e., 

^max = max { > R V } < $ . (4-12) 


where R $ is the sum of the normalized absolute residuals across all the computational nodes, 
namely 


R<v 


I 


£(A,<P,) + SAV-A,<P p 

J = E,W.N,S 


F.0.. 


(4-13) 


and 0 = 1, U and V for the mass equation, momentum equations in axial and radial directions, 
respectively. F in is the mass-flow rate at the inlet, and, O in takes the inlet axial mean velocity for 

the momentum equations. Note that the azimuthal velocity W does not appear in the mass 
conservation equation for axisymmetric flows. Therefore, it is not used in the criteria for 
convergence (4-12). 
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4.2.3 Cases studied Swirling flows with a reverse flow region provide further challenges for 
turbulence modeling. To rigorously test our turbulent models, the experiments of Weske and 
Sturov (1974) and Kitoh (1991) are selected for model validations. The reverse flow region 
observed in Kitoh's (1991) experiments starts near the inlet and has a length about 40 diameters of 
the pipe. According to the definition of swirl number (4-1), the reference or entrance swirl number 
S 0 is about 1.3 in Weske and Sturov's case and approximately 1 in Kitoh's (1991) case; the 

Reynolds number (based on the pipe diameter and averaged mean axial velocity) is about 30,000 
and 50,000, respectively. As mentioned before, due to different ways of generating swirl, the 
mean azimuthal velocity at the inlet of the pipe is forced-vortex dominated in the Weske and 
Sturov's case, whereas it is free-vortex-type in the Kitoh's case. As a result, the azimuthal 
velocity profiles evolve differently downstream. In Weske and Sturov's case, the dominated 
forced-vortex region reduces and the free-vortex region grows as swirl decays, with maximum 
azimuthal velocity shifting to the centerline of the pipe. In Kitoh's case, the free-vortex region is 
dominated in the reverse flow region, and initial swirl decay is only associated with free vortex 
motion. 

4.2.4 Comparisons with data The computational domain is one radius in the radial direction 
and 100 diameters in the axial direction. Two sets of grid points, 102x86 and 51 x 56 in the 
axial and radial directions (Figs. 4.4a - 4.4b), are tested. The first grid shown in Fig 4.4a is 
clustered in the inlet region, near-wall region and core (near centerline) region. The second grid 
shown in Fig. 4.4b is also clustered in near-wall region, but less dense in the inlet and core 
regions. For the NNWRS and NWRS models, both grids give almost the same predictions of the 
mean and turbulence flow fields. Thus, only the results from the grid 51 x 56 are presented here. 

The numerical results are presented first for the Weske and Sturov's (1974) case and then for 
the Kitoh's (1991) case. Fig. 4.5 shows the streamline defined on the (x, r)-plane given by the 
NNWRS model. Similar flow pattern given by the NWRS model (no shown) is observed. The 
predicted flow pattern indicates no reserve flow occur in the flow region, which agrees with the 
experimental observations. Owing to small radial velocity V, the streamlines on ( x , r)-plane are 
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almost parallel to the centerline. The predicted swirl intensity along the pipe is compared with the 
data in Fig. 4.6. Both models give about the same swirl decay rate, which gradually decreases 
along the pipe and approaches to a constant in the downstream. The predicted swirl number agrees 
very well with the data up to x/D = 20, but it is slightly overestimated at x/D = 50. Further 
downstream, no data are available for comparisons. 

Figs. 4.7-4.9 compare model predictions of the mean flow and turbulent fields with available 
data at measured locations x / D = 5. 1 , 20 and 50 (note the first measured location, corresponding 
to x -0 in the computational domain, is x / D = 0.35). Both models give almost the same 
results. The agreement between the predicted mean axial velocity U and the data is very good 
(Figs. 4.7a, 4.8a, and 4.9a). The mean axial velocity in the core region increases with increasing 
axial distance (x), and its profile becomes almost uniform across the entire cross-section (except 
the near-wall region) in the downstream. The predicted azimuthal velocity W also agree well with 
the data (Figs. 4.7b, 4.8b, and 4.9b), although it gives a slightly higher vorticity strength in the 
forced-vortex region and its peak location is shifted a little bit toward to the centerline. The peak 
value of W decreases due to the decay of swirl along the pipe, and W in the downstream becomes 
more uniform (cf. Fig. 4.7b with 4.9b). Both models essentially repeat the experimental findings 
that the forced-vortex region shrinks and the vorticity strength weakens as the distance from the 
inlet increases, whereas the free-vortex region increases with increasing distance. Predicted flow 
skewness near the wall at different measured locations is shown in Fig. 4.10 by plotting U as a 
function of W . Very close to the wall, the ratio U /W is a constant, which becomes larger as x 
increases. 

The agreement between the predicted turbulent quantities (kinetic energy k and root-mean- 
squared turbulent fluctuations u rms , v rns , and w rms ) and the data is not very satisfactory (Figs. 
4.7c-f and 4.9c-f; note that no experimental data are available at x/D = 20 for k, u rms , v rms , and 
w^). Near the inlet, at x / D = 5.1, the models predict the trends of these turbulent quantities 
reasonably well, and the agreement between the model predictions and the data in the core region is 
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acceptable (Figs. 4.7c-f). Further downstream, at x/D = 50, the model predictions agree well 
with the data in the region covering half the radius from the wall, but are underestimated in the core 
region (Figs. 4.9c-f). 

For the Kitoh’s case, the models are able to capture the recirculation in the core region near 
the inlet reported by Kitoh (1991) (Fig. 4.1 1). However, the length of the reverse flow region 
predicted by the models is only about 7 diameters, which is much smaller than the observed length. 
This accounts for the major discrepancies between the predicted and measured flow fields 
discussed below. 

The predicted swirl number along the pipe is compared with the data and the empirical 
formula given by Kitoh (1991) (Fig. 4.12). The agreement among them is excellent. Swirl 
attenuates exponentially along the pipe with a constant decay rate. Figs. 4. 1 3 - 4. 1 7 compare the 
model results with the experimental data at following measured locations, x! D= 12.3, 19.0, 
25.4, 32.4 and 39.0 (note that the first measured location, corresponding x = 0 in the 
computational domain, is x / D = 5.7). At each location, six plots are presented for the mean 
velocities U and W , turbulence kinetic energy k the root-mean-squares of turbulent fluctuations 
u r ms’ v rms > ar >d w ms ■ Overall, both models give almost the same predictions for mean flow and 
turbulent field (Figs. 4.13-4.17). Reverse axial velocity (U < 0) (predicted by the models) is 
observed only at x/ D= 12.3, the closest location to the inlet (cf. lines in Fig. 4.13a with 4.14a- 
4.17a). This is consistent with the flow pattern shown in Fig. 4.11. At each location, the 
predicted mean velocities follow the data well in the near-wall region, whereas they depart from the 
data in the core region of the pipe with U being overpredicted and W being underpredicted. The 
disagreement becomes worse as the location is further downstream. Fig. 4.18 shows the predicted 
flow skewness in the near wall region. Again, very close to the wall, U is a linear function of W 
with a slope increases as the distance from the inlet increases. 

For turbulence kinetic energy k and root-mean-squared turbulent fluctuations u rms , v rm( , and 
w rms , g°°d agreement between the model predictions and the data is found in the near-wall region, 
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while discrepancy is observed in the core region of the pipe (Figs. 4. 13c-f, 4.14c-f, 4.15c-f, 
4.16c-f; 4.17c-f). Both models perform consistently in predicting the mean velocities and in 
predicting the turbulence variables. Note that in this case the turbulent stresses are quite irregular 
in the core region due to the recirculation. This case is a real challenge for any turbulence model. 

4.3 Turbulent Developing Flow in a Rotating Pipe 

4.3.1 Background Flows in a rotating pipe can be classified as a subset of swirling flow 
discussed in the previous section because they are also subject to extra shear strain associated with 
the mean azimuthal velocity component. However, swirl intensity does not decay in rotating pipe 
flow because angular momentum is continuously imparted into the flow. In the downstream, the 
flow approaches to an equilibrium state as the profile of the mean azimuthal velocity approaches to 
a rigid-body-rotation distribution. 

The rotation effect on a fully-developed pipe flow was studied by White (1964) and 
Murakami and Kikuyama (1980). Both studies showed that rotation stabilizes the turbulent flow 
field and reduces the flow resistance with a rate increasing as the rotation increases. In other 
words, rotation promotes the relaminization of a fully-developed turbulent flow. Experimental 
study of the rotation effect on developing turbulent flow in a rotating pipe was carried out by 
Kikuyama et al. (1983). They found out that when flow enters a rotating pipe with a uniform 
entrance velocity, it is affected by both destabilizing effect created by a large circumferential shear 
strain and stabilizing effect due to centrifugal force. Near the entrance, the wall boundary layer is 
very thin, and the mean azimuthal velocity has to decrease rapidly from the circumferential velocity 
of the wall to zero outside the boundary layer. The flow near the wall is subject to a very high 
mean circumferential shear strain, which results in great enhancement of the turbulence production. 
Therefore, the near-wall flow in the upstream region is destabilized by the dominated destabilizing 
effect. As flow moves to the downstream, the boundary layer becomes thicker and thicker, and the 
intensity of the destabilizing effect decreases due to the decrease of the circumferential shear strain. 
Far downstream, the flow becomes fully developed, and the fluid inside the rotating pipe rotates as 
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a rigid body. Since the rigid-body motion does not give rise to additional turbulence production, 
the flow is stabilized by the stabilizing effect of centrifugal force. It follows that as the flow 
develops along the rotating pipe, the destabilizing effect gives away to the stabilizing effect, and the 
flow will pass through a region where both effects are equally important. Hence, the flow in a 
rotating pipe is very complicated and provides a severe test for turbulent modeling. 

Flow characteristics associated with developing turbulent flows in a rotating pipe are 
commonly found in the inlet part of fluid machines, heat exchangers, and cooling system of the 
rotors. Therefore, correct modeling of rotating turbulent flows is very important for the design and 
development of any rotating machines. Attempts have been made to model developing turbulent 
flow in a rotating pipe through the streamline curvature modification with partial success. Most 
recently Yoo et al. (1991) used a near-wall Reynolds stress model based on the LRR pressure- 
strain model discussed in Section 2 to model the rotating pipe flow. Their model gave fairly good 
predictions in comparison with Kikuyama et al.'s (1983) data. In this section, the performance of 
the NNWRS model will be assessed by comparing the model predictions with Kikuyama et al.'s 
(1983) data and those given by the NWRS model and Yoo et al.'s model. 

4.3.2 Numerical considerations The governing equations for rotating pipe flow are the same 
as those for swirling flow. At the inlet of the flow domain, the mean axial velocity is determined 
from the measurement, and the turbulence quantities are estimated from the inlet mean axial 
velocity U in . More specifically, the turbulent kinetic energy k in and normal stresses 
uuin, vv,>j,and wwi n are specified according to the following formula: 


k - IT ! 2 
K in lu in ’ 

(4- 14a) 

UUin - k in , 

(4- 14b) 

VVin=0.6k in , 

(4- 14c) 

WWin=0Ak in , 

(4- 1 4d) 
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where / is the measured turbulence intensity at the inlet, the dissipation rate is estimated according 
to the formula (4-3), and the turbulent shear stresses are approximated by Boussinesq 
approximation. Negligible diffusion conditions in the axial direction are imposed at the outlet of 
the computational domain, and axisymmetry condition (4-4) is used at the centerline. No-slip 
condition requires that at the moving wall W equal to the circumferential velocity of the pipe 
surface W 0 , and other variables, except mean pressure and dissipation rate, be zero. The 

dissipation rate £ at the moving wall is given by (4-5). Same numerical procedure discussed in 
section 4.2 is adopted here for the calculation of developing flow in a rotating pipe. 

4.3.3 Results and discussion The developing pipe flow experiments of Kikuyama et al. 
(1983) are used to validate our models. Three sets of experimental data with rotation number 
N = W 0 / U m = 0, 0.5 and 0.83 are available, where U m is the averaged mean axial velocity at the 

inlet. The corresponding swirl number S is 0, 0.25, and 0.415, respectively. The entrance 
velocity is almost uniform with Reynolds number Re = U m DI u = 6x 10 4 . The turbulence 

intensity I at the entrance was about 0.3%. Both mean and turbulence flow quantities were 
measured at several locations. 

The computational domain again is 0.5D (r)x 100D (x) with a 51 x 56 grid in radial and 
axial directions. The grid is found to be able to give grid independent results. Numerical 
computations are carried out for three cases with rotation numbers N = 0, 0.5 and 0.83 in 
Kikuyama's (1983) experiments. The results for the most severe case, namely N = 0.83 case is 
presented below to show the models' ability to predict developing turbulent flow in a fast rotating 
pipe. For other two smaller rotation number cases, model predictions are in better agreement with 
the measurement. The rotation effects on the turbulent kinetic energy is shown as N increases 
from 0 to 0.83. 

Model predictions of the mean axial and azimuthal velocity profiles are compared with the 
data at the first and last measured locations x/ D = 2.7 and 28.5 in Fig. 4.19a-4.19b. All three 
models give good predictions of the mean velocity profiles, and the results given by the NNWRS 
and NWRS models are very close. For the axial velocity U , the NNWRS and NWRS models 
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give better predictions at the first measured location x / D = 2.7 than at the last measured location 
x / D - 28. 5 . For the azimuthal velocity W , the agreement between the predictions given by the 
NNWRS and NWRS models and the data is almost the same at these two locations. Our model 
predictions of U and W follow the trend of the data with a small deviation. Yoo et al.'s (1991) 
model gives slightly better predictions of U and W in terms of magnitude, but not slope. The 
measurements in Fig. 4. 19 indicate that the boundary layer becomes more developed downstream, 
and the mean axial velocity profile at x / D = 28.5 is close to that in a fully-developed pipe flow. 
These features are essentially captured by all three models. 

The comparisons for Reynolds stress components at x / D = 2. 7 and x ID = 28.5 are 
presented in Figs. 4.20-4.25. Turbulent fluctuations (or normal stresses) and shear stresses 
spread out into the core region as the destabilizing effect associated with the circumferential shear 
strain decreases and the stabilizing effect due to centrifugal force becomes more dominated (cf. 
panel a to panel b in Figs. 4.20-4.25). In the near-wall region, all the models overpredict the 
turbulent fluctuations (normal stresses) u rms v rms and w rms , except u rms at x! D = 2.7 where 
the peak value is predicted well by the models (Figs. 4.20 - 4.22). Yoo et al.'s model much more 
overpredicts the peak values of the v rms than the NNWRS and NWRS models do (see Fig. 4.21a 
and 4.21b). Overall, the predicted turbulent fluctuations by the NNWRS and NWRS models agree 
slightly better with the measurements than those by Yoo et al.'s model do. The components of the 
shear stresses are generally overpredicted by the NNWRS and NWRS models (Figs. 4.23 - 4.25) 
except for vw at x/ D = 28.5 and uw at x I D = 2.7 (Figs 4.24b and 4.25a). The 
performance of Yoo et al.'s model varies depending on the individual component and the measured 
location. It gives good predictions of each shear strain component at one location, but either 
overpredicts or underpredicts the same component at the other location. 

To show how the mean velocity profiles evolve along the pipe from the upstream to the 
downstream, we plot the axial and azimuthal velocity components at four different measured 
locations: x / D = 2.7, 9.7, 15.5 and 28.5 in Figs. 4.26a - 4.26b. Note that only the results of 
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the NNWRS model and the measurements are compared. Both the axial and azimuthal velocities 
have similar profiles at these four different locations. As the flow moves to the downstream, the 
boundary layer thickness increases and the variations of U and W across the boundary layer 
become less rapid. The NNWRS model follows the trends of the measurements quite well but 
with a small offset. Distributions of the turbulent kinetic energy k at the same locations are shown 
in Fig. 4.26c. As expected, the turbulent kinetic energy k spreads out into the core region 
resulting from the increase of the stabilizing effect and the decrease of the destabilizing effect in the 
downstream. The NNWRS model replicates this observed phenomenon, but it overestimates the 
maximum kinetic energy and underestimates the spreading rate. 

The effects of rotation number N on the mean axial and azimuthal velocity profiles at the last 
measured location x! D — 28.5 are shown in Figs. 4.27a - 4.27b. The NNWRS model gives 
better predictions for the moderate rotation number N = 0.5 case than for the high rotation number 
N = 0.83 case. Both the measurements and model predictions show that the mean axial and 
azimuthal velocity components far away from the inlet are not very sensitive to the rotation rate. 
Finally, the effects of rotation number N on evolution of the turbulent kinetic energy k are shown 
in Fig. 4.28. At the upstream location x / D = 2.7, the NNWRS model shows an increase of k 
with the rotation rate. This is because more turbulent production associated with the more intense 
shear strain is generated as rotating rates increases. At x / D = 15.5, the total kinetic energy for 
different rotation numbers are close to each other since the stabilizing and destabilizing effects are 
well balanced there. The drop in k at the downstream location xl D = 28.5 indicates the 
stabilizing effect there is dominated, resulting in the decrease of the turbulent kinetic energy with 
increasing rotation rate. These predictions are consistent with the experimental measurements 
given by Kikuyama et al. (1983). 

4.4 Fully-Developed Turbulent Flow in a Rotating Channel 

4.4.1 Background For fully developed flow in a spanwise rotating channel with constant 
angular velocity £2 (Fig. 4.29), the system rotation gives rise to two additional body forces: 
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Coriolis force and centrifugal force. However, only the Coriolis force is dynamically important 
because the centrifugal force term can be absorbed in the effective pressure 


p eff = p--p(nr y 


(4-15) 


where P is the stationary pressure and r denotes the distance from the rotating axis. The strength 
of the Coriolis force can be measured by a rotation number 


Ro = 


-® , or = ® 

U m 


(4-16) 


which is the relative strength of the Coriolis force to the inertial force. Here U m is the bulk mean 
velocity, h is the half the channel width, and u T is the global friction velocity. 


The flow under consideration is driven by an imposed mean pressure gradient, 

iP eff _ P»t 
dx h 


(4-17) 


in the streamwise (x) direction. Thus, the turbulence can be treated as homogeneous in the 
streamwise (x) and spanwise (z) directions. In other words, the mean properties and turbulence 
statistics vary only in the transverse (y) direction. With the use of effective pressure, the 
governing equation for the mean flow is identical to that for non-rotating channel flow, i.e., 


1 dp eff d( dU — 

— = — v uv 

p dx dy \ dy 


(4-18) 


It follows from (4-17) and (4-18) that 

u l=l( u l + u l>) > 


(4-19) 


where u Ts and u ^ denote the local friction velocity at the suction side and pressure side, 

respectively. When Q>0 (<0), y/ h = +\ is the suction (pressure) side and y/h = - 1 is the 

pressure (suction) side. Without system rotation, the mean flow is symmetric with respect to the 
(x, z)-plane at y = 0 and u r = u Ts = u rp . 
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In the transport equations for the individual Reynolds stress components, the production 
terms associated with mean shear and rotational stress generation are 




+ UjU k 


dUj 

dx k 


(4-20a) 


Gy - ~2Q(. 


\u j u m e ikm + Lt t li rn t ’jkin ] 


(4-20b) 


respectively. The components of and Gy are given in Table 4.3. Since G kk vanishes, no 
additional turbulence energy is generated directly through the rotational production term Gy. 

Pioneering work on fully developed flow in a channel rotating about a spanwise axis was 
carried out by Johnston et al. (1972). Their measurements covered a wide range of rotation 
numbers, Ro, from 0 up to 0.21 with Reynolds number Re ranging from 11500 to 35000. Here 
the Reynolds number is defined based on the bulk mean velocity U m and the channel width 2h as 
Re = 2 U m h / v. The Coriolis force was found to affect both local and global stability of the flow. 
Three stability-related phenomena caused by the Coriolis force were observed or inferred: (1) it can 
change the streak bursting rate in the wall-layers; (2) it can suppress the turbulence production near 
the suction (stabilized) side; (3) it can develop a large-scale roll cell on the pressure (destabilized) 
side. 

Theoretical and experimental studies were performed by Nakabayashi and Kitoh (1996) on 
low-Reynolds-number, fully developed turbulent flow in a rotating channel. Low-Reynolds- 
number flow is found to be more strongly affected by the Coriolis force than high-Reynolds- 
number flow. The Coriolis force affects the logarithmic layer and the core region. By dimensional 
analysis, they deduced that Reynolds number R e T = u r h/ v and Ro are the two parameters 
determining the overall flow structure. The ranges of these two parameters in the experiments they 
investigated are 28 < Re T < 1 55 ( 1700 < Re < 10000) and 0 < Ro < 0.055. 

Laboratory investigations of turbulence in rotating reference systems are more difficult to 
accomplish than most other experiments, simply because the flow apparatus has to be mounted on 
a rotating turntable. While fully developed flows are relatively difficult to realize in the laboratory, 
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especially in a rotating system, such flows have been attractive objects for numerical study in the 
past two decades. Kim (1983) was the first to perform large-eddy simulations (LES) for rotating 
channel flow with a moderate rotation number Ro = 0.068, which revealed how the Coriolis force 
introduced by the system rotation can stabilize or destabilize turbulence on different sides of the 
plates. Miyake and Kajishima (1986) also carried out large -eddy simulations for rotation number 
up to 0.2. They presented various statistical turbulence quantities and concluded that in the near- 
wall region, the Coriolis force enhances sweep and ejection on the pressure side, while reduces 
them on the suction side. However, the rather course grid used made their results less reliable and 
the conclusions less convincing. 

Kristoffersen and Anderson (1993) performed direct numerical simulations (DNS) of fully- 
developed rotating channel flow for rotation number Ro up to 0.50 at Re = 5800. Their 
simulation results showed that with increasing rotation the damping and augmentation of 
turbulence along the suction side and pressure side, respectively, become more significant, 
resulting in highly asymmetric profiles of mean velocity and turbulent Reynolds stresses. The 
mean velocity profile exhibits an appreciable region with slope 2Q, in accordance with the 
experimental observations of Johnston et al. (1972). At Ro = 0.50 the Reynolds stresses vanishes 
in the vicinity of the suction (stabilized) side. Because the gird they used, 1 28 x 1 28 x 1 28 , is find 
enough to resolve all essential scales of the low-Reynolds-number turbulence, their rather complete 
data become very useful to verify turbulent models for rotating turbulent flow. 

A Reynolds stress model is required to model rotating channel flow because of high 
anisotropic turbulence caused by the Coriolis force. Second-order turbulence models have been 
used to predict rotating channel flow recently. Launder et al. (1987) and Shima (1993) obtained 
good agreement with the experimental results of Johnston et al. (1972) and large-eddy simulations 
of Kim (1983). Tselepidakis (1991) made an attempt to model the fully developed rotating channel 
flow studied by Kristoffersen and Anderson (1993). Very good agreement was obtained with low 
rotation number Ro = 0.05. However, no results were reported for higher rotation number due to 
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numerical stability. He suggested that the problem could be overcome by reconstructing a new 
dissipation rate tensor £ y . 

4.4.2 Results and discussion In the present study, the cases studied by Kim (1983) and 
Kristoffersen and Anderson (1993) are selected to test our near wall hnodels' ability to predict fully 
developed channel flow with spanwise rotation. In Kim's case, the rotation number is Ro = 0.068 
(Ro t = 1.47) and the Reynolds number is Re c =2U c h/ v = 13800 (where U c is the centerline 
mean velocity) (Re T = 640); in Kristoffersen and Anderson's case the rotation number Ro( Ro r ) 
varies from 0 to 0.50 (7.55) at Reynolds number Re = 5800 (Re T = 194). The iterative scheme 
for the fully developed channel flow in Section 3 is adopted here for the calculations of the fully 
developed flow in a rotating channel. 

The numerical results are presented first for Kim's (1983) case, and then for Kristoffersen 
and Anderson's (1993) case. Note that the channel in these two cases rotates in opposite 
directions, i.e. Q>0 in Kim's case, whereas £2<0 in Kristoffersen and Anderson's case. 
Accordingly, y/h = - 1 is the pressure side and y/h = + 1 is the suction side in Kim's case, 
whereas y/h = - 1 is the suction side and y/h = +\ is the pressure side in Kristoffersen and 
Anderson's case. The global wall friction velocity u r is used to normalize the flow variables when 
variations across the entire channel -\<y / h<+l are presented, whereas the local friction 
velocities u ^ and u Xs are used in the scaling of computed results whenever the inner coordinate 
y + near the pressure and suction sides, respectively, labels the abscissa. In the computations, two 
damping functions, (2-29a) for Kim’s case and (2-29b) for Kristoffersen and Anderson's case, are 
used in the NWRS model due to the significant difference in Reynolds numbers, whereas same 
damping function (2-33) is used in the NNWRS model for both cases. 

Figs. 4.30a - 4.30b compare the predicted mean velocity with the LES data given by Kim 
(1983) in wall coordinates. The LES data show that the normalized mean velocity U + lies above 
and below the semi-logarithmic law U + = 2.51ny + +5.5 in the region near the centerline on the 
suction side (Fig. 4.30a) and pressure side (Fig. 4.30b), respectively. On the suction side 
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( y / h = +1), the wake component is increased due to the decrease of local friction velocity u TS (i.e. 
u rs <u r ). On the pressure side (y/h = -\), the increased local friction velocity (i.e. u ^ > u T ) 
suppresses the wake component in the velocity profile. Both the NNWRS and NWRS models are 
basically able to capture these features; the NWRS model is in better agreement with the LES data 
than the NNWRS model does. 

Fig. 4.31a plots the normalized mean velocity U / U c distribution between the plates. The 
model results are essentially identical and follow the LES data quite well. The data and the model 
results show a slight asymmetry in the mean velocity profiles due to the presence of rotation. Both 
models give about the same shear stress uv predictions as shown in Fig. 4.31b. Fig. 4.31c shows 
that the peak value decreases on the suction side (y/h = +1) and increases on the pressure side 
(y/h = - 1) in comparison with the non-rotating case (no shown). Note that although the 
rotational production term G \j does not have a direct contribution to the generation of turbulent 
kinetic energy, it affects the mean shear production term through changing the turbulent stress 
distributions and the mean velocity gradient, resulting in the redistribution of the kinetic energy k + 
between the plates. The NWRS model predicts a much higher peak value of k + on both sides 
compared with the NNWRS model. Away from the walls, both models give almost the same 
predictions, which are slightly overestimated on the suction side and underestimated on the 
pressure side. In Fig. 4.3 Id, the LES shows that in the core region, turbulence is more isotropic 
on the pressure side than on the suction side. Both models also show that turbulent fluctuations 
are more isotropic (i.e. u+ ms , v+ ms and are closer to each other) on the pressure side. 

However, except for u+ ms on the pressure side, model predictions of turbulent fluctuations do not 
agree well with the LES data. Overall, the models seem to have the ability to replicate the rotation 
effect, but they are less responsive to the influence of the rotation. Further verification for the 
models is needed since in this case the rotation number Ro = 0.068 is relatively small compared 
with its high Reynolds number Re c = 13800. 
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According to Nakabayashi and Kitoh (1996), with the same rotation number, the rotation 
effect is more influential in low-Reynolds-number flow than in high-Reynolds-number flow. To 
further examine the rotation effect in the model, computations are carried out for the low-Reynolds- 
number case of Kristoffersen and Anderson with rotation numbers Ro = 0, 0.01, 0.05, 0.10, 
0.15, 0.20, and 0.50 (these values correspond to Ro r = 0, 0.15, 0.755, 1.51, 2.265, 3.02, and 
7.55, respectively). The calculation for rotation number Ro = 0 is to test the code and provide a 
reference to identify the rotation effects (Fig. 4.32-4.33). The NWRS and NNWRS models repeat 
their performance in fully developed channel flow without rotation. Note that all flow variable 
profiles are symmetric about the center plane y = 0 (Fig 4.33). In the following presentation, the 
mean velocity profiles are shown in the wall and global coordinates for each rotation number, 
whereas the turbulence quantities are plotted only in the global coordinate. 

Figs. 4.34a - 4.34b plot the mean velocity profiles in the wall coordinates for Ro = 0.01. 
According to the DNS, the effect of the Coriolis force at this low rotation number can be 
considered as a small perturbation to the non-rotating case, and the velocity profiles in the core 
region on different sides deviate slightly from the log law U + = 2. 5 In y + + 5.5 . The deviation on 
the pressure side is more noticeable than that on the suction side. Our model predictions agree very 
well with the DNS data on the suction side (Fig. 4.34a). On the pressure side, the NNWRS model 
predicts a slight larger mean velocity in the core region compared with the DNS data, whereas the 
predictions given by the NWRS model are closer to the DNS data. Predicted mean velocity and 
turbulence quantities are further compared with the DNS data in the global coordinate (Fig. 4.35). 
The agreement between our model predictions and the DNS data is very close to that in the non- 
rotating case (cf. Fig. 4.35 with its counterpart Fig. 4.33). The weak rotation in this case only 
causes a small perturbation from the non-rotating state. 

The predicted mean velocity profiles and the DNS data for Ro = 0.05 are shown in Fig. 
4.36. The DNS data show that the rotation effect becomes more pronounced, namely the wake 
component is appreciably enhanced on the suction side and is eliminated on the pressure side. This 
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is accompanied by the rather significant decrease and increase of the wall fiction velocity on the 
suction and pressure sides respectively. The velocity profiles predicted by the NNWRS and 
NWRS models, however, only slightly deviate from the law of the wall. Compared to their 
predictions for /to = 0.01 (cf. Fig. 4.36 with Fig 4.34), the models are not sensitive to the 
increase of rotation. This insensitivity is also observed when mean velocity and turbulence 
quantities are plotted in the global coordinate (compare the model predictions in Fig. 4.37 with Fig. 
4.35). The DNS data show that the rotation already has a significant influence on the turbulent 
field at this rotation number. The peak of the shear stress uv increases on the pressure side but 
decreases on the suction side (cf. Fig. 4.37b with Figs. 4.33b and 4.35b). This is because the 
rotational production term G ( y (see Table 4.3) in the wv-transport equation is positive near the wall 
and therefore tends to reduce the negative value of uv on the suction side and increase it on the 
pressure side. Away from the walls, the viscous effect is negligible, and the shear stress profile 
becomes linear (see (4-18)). The peak value of the turbulent kinetic energy drops on the suction 
side (y / h = -1) and rises on the pressure side (y / h = + 1) quite significantly (Fig. 4.37c); the 
peaks of turbulent fluctuations v+ ms and w^ ms in the vicinity of the suction side tends to disappear 
(Fig. 4.37d). The models do not predict such features. The rotation number Ro = 0.05 in this 
case is close to Ro = 0.068 in Kim's case. Comparing the LES data in Figs. 4.30-4.31 and the 
corresponding DNS data in Figs. 4.36 - 4.37, we find that the rotation effect is more influential in 
this low-Reynolds-number case (Re T = 194) than in the high-Reynolds-number Kim's case 
(Re T = 640), which is in accordance with Nakabayashi and Kitoh (1996). 

As the rotation number increases up to Ro = 0.20, comparisons between the models’ results 
and the data are qualitatively similar to those for Ro = 0.05 with some quantitative differences. 
The comparisons for mean velocity profiles in the wall coordinates are presented in Figs. 4.38, 
4.40, and 4.42 for rotation number /to = 0.10, 0.15 and 0.20, respectively, whereas those for 
turbulence quantities can be found in Figs. 4.39, 4.41, and 4.43. 
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At Ro = 0.50, the difference in the mean flow and turbulent field on the suction side and 
pressure side becomes striking. The DNS shows that the mean velocity profile is close to the 
linear law U + - y + on the suction side, whereas an anomalous behavior is observed on the 
pressure side: the profile diverges from the linear law at y + ~5. The NNWRS and NWRS 
models seem to be able to capture these features, although the predicted profiles deviate from the 
data. Unlike those lower rotation number cases (Ro <0.5) (see Figs. 4.34b, 4.36b, 4.38b, 
4.40b, and 4.42b), the difference in the predicted mean velocity profiles is quite large on the 
pressure side in this case. The NWRS model results are in a better agreement with the data. In the 
global coordinate (Fig. 4.45a), a linear region in the DNS mean velocity profile is observed. Both 
models are able to predict this characteristic flow region, though the predicted slope is smaller than 
the expected 20 as indicated by the data. The predicted shear stress uv shows a better agreement 
with the data in this case than in all lower rotation number cases except for Ro = 0.01 case (cf. 
Fig. 4.45b with Figs. 4.37b, 4.39b, 4.41b, and 4.43b). More pronounced rotation effects are 
also observed in the kinetic energy and turbulence fluctuations at Ro = 0.50 (Figs. 4.45c - 4.45d). 
According to the data, on the suction side, the strong rotation effect has eliminated the peaks of the 
kinetic energy k + and the normal stress in the wall normal direction , and tends to eliminate 
the peaks in the other two normal stresses u^ ms and w+ ms . In the core region, is dominant 
among the three normal stress components. The NNWRS and NWRS models have shown a more 
significant rotation effect at Ro = 0.50 in comparison with the lower Ro cases. More asymmetry 
due to rotation in the predictions is observed at this rotation number. However, the rotation effect 
predicted by the models are still much weaker than that shown in the DNS. For instance, the peaks 
in the kinetic energy k + and normal stresses u+ ms , v+ ms , and w+ ms near the suction side still 
remain. 

To show more clearly the variations of the predicted mean velocity as the rotation rate 
increases, we plot the mean velocity profiles predicted by the NNWRS model for different Ro in 
the wall coordinates together in Fig. 4.46a and Fig. 4.46b. The mean velocity profiles on the 
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suction side lie above the profile for the non-rotating case in the semi-logarithmic region. The 
profile for small Ro still exhibits a logarithmic region. As Ro increases, the profile has a more 
laminar-like shape and approaches the linear law U + = y + . Correspondingly, the mean velocity 
profiles on the pressure side lie below the profile for Ro = 0. These profiles still exhibit the 
characteristic turbulent semi-logarithmic shape, but with a slope that decreases as Ro increases. At 
the highest rotation rate ( Ro = 0.50), an anomalous behavior is observed. 

Finally, Fig. 4.47 shows the effect of rotation on local wall friction velocities u Ts and u Tp . 
The NNWRS and NWRS models follow the trend indicated by the DNS data, namely the wall 
friction velocity decreases on the suction side and increases on the pressure side as the rotation rate 
increases. The model results, however, show a weaker rotation effect on the wall friction 
velocities. 

4.5 Concluding Remarks 

The performance of the NNWRS model in predicting complex turbulent flow with swirl and 
rotation is examined in this section. More specifically, we apply the model to swirling pipe flow, 
developing rotating pipe flow, and fully developed channel flow with spanwise rotation and 
compare the model results with available experimental and LES/DNS data. The complexity of 
these flows provides severe tests for our models. In swirling pipe flow, the NNWRS model gives 
reasonable predictions for flow with no reverse flow region, whereas for flow with a reverse flow 
region, the model predicts a much smaller reverse flow region than that observed in the 
experiments. In developing flow along a pipe rotating around its axis, the NNWRS model is able 
to capture the flow pattern subject to both destabilizing effect resulting from the mean 
circumferential shear strain and stabilizing effect due to the centrifugal force. In the fully 
developed channel flow with spanwise rotation, at low rotation rate, the model predictions agree 
with the LES/DNS data. As the rotation rate increases, the model underestimates the rotation 
effect, which causes the augmentation and damping of the turbulence along the pressure and 
suction sides, respectively. The NWRS model behaves similarly as the NNWRS model in all the 
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cases considered. This implies that the near-wall modifications in the turbulence models become 
less important for complex flows such as swirling and rotating flows. Further improvement on the 
modeling is needed for a better agreement with the experimental data and the direct numerical/large- 
eddy simulations, in particular in the reverse flow region of the swirling pipe flow and semi- 
logarithmic region of the rotating channel flow where the Coriolis force is important. In the next 
section, the NNWRS model is further examined for turbulent flows with complex geometry. 
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Table 4. 1 Comparison of shear strain rates in swirling and parallel flows 
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Table 4.2 Swirling flow characteristics 
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Table 4.3 Production terms in fully developed channel with spanwise rotation 
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Fig. 4.1 Typical azimuthal velocity profile in a swirling pipe flow. 


Control volume 



Fig. 4.2 Control volume centered at grid point P . 
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Control volume for U 



S 


Fig. 4.3a Control volume for mean velocity in the axial direction ( U). 
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Control volume for V 


Fig. 4.3b Control volume for mean velocity in the radial direction ( V). 
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Fig. 4.6 Comparison between the predicted and the measured swirl number for 

Weske and Sturov's (1974) case. 
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Fig. 4.7a Comparison between the predicted and measured mean axial velocity profiles at 

measured location x / D = 5. 1 . 



Fig. 4.7b Comparison between the predicted and measured mean azimuthal velocity profiles at 

measured location x / D = 5.1. 
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Fig. 4.7c Comparison between the predicted and measured turbulent kinetic energy at measured 

location x / D = 5. 1 . 



Fig. 4.7d Comparison between the predicted and measured turbulent fluctuation in the axial 

direction at measured location x! D — 5.1 . 
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Fig. 4.7e Comparison between the predicted and measured turbulent fluctuation in the radial 

direction at measured location x / D = 5. 1 . 



Fig. 4.7f Comparison between the predicted and measured turbulent fluctuation in azimuthal 

direction at measured location x / D = 5. 1 . 
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Fig. 4.8a Comparison between the predicted and measured mean axial velocity profiles at 

measured location x / D = 20.0 . 



Fig. 4.8b Comparison between the predicted and measured mean azimuthal velocity profiles at 

measured location x / D = 20.0. 
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Fig. 4.9a Comparison between the predicted and measured mean axial velocity profiles at 

measured location x / D - 50.0. 



Fig. 4.9b Comparison between the predicted and measured mean azimuthal velocity profiles at 

measured location x 7 D - 50.0. 
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Fig. 4.9c Comparison between the predicted and measured turbulent kinetic energy at measured 

location x I D = 50.0 . 



Fig. 4.9d Comparison between the predicted and measured turbulent fluctuation in the axial 

direction at measured location x ID — 50.0. 
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Fig. 4.9e Comparison between the predicted and measured turbulent fluctuation in the radial 

direction at measured location x / D = 50.0. 



Fig. 4.9f Comparison between the predicted and measured turbulent fluctuation in azimuthal 

direction at measured location x / D- 50.0 . 
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Fig. 4.12 Predicted swirl number for Kitoh’s (1991) case and its comparison with the experimental 
data and the empirical formula given by Kitoh (1991). 
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Fig. 4.13a Comparison between the predicted and measured mean axial velocity profiles at 

measured location x / D = 12.3. 



Fig. 4.13b Comparison between the predicted and measured mean azimuthal velocity profiles at 

measured location at x / D = 12.3. 
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Fig. 4.13c Comparison between the predicted and measured turbulent kinetic energy at measured 

location jc / D = 12.3. 



Fig. 4. 13d Comparison between the predicted and measured turbulent fluctuation in the axial 

direction at measured location x / D = 12.3. 
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Fig. 4.13e Comparison between the predicted and measured turbulent fluctuation in the radial 

direction at measured location x / D= 12.3. 



Fig. 4. 1 3f Comparison between the predicted and measured turbulent fluctuation in azimuthal 

direction at measured location x / D = 1 2.3. 
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Fig. 4.14a Comparison between the predicted and measured mean axial velocity profiles at 

measured location x / D = 19.0 . 



Fig. 4.14b Comparison between the predicted and measured mean azimuthal velocity profiles at 

measured location x / D — 19.0 . 
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Fig. 4.14c Comparison between the predicted and measured turbulent kinetic energy at measured 

location x / D= 19.0. 



Fig. 4.14dComparison between the predicted and measured turbulent fluctuation in the axial 
direction at measured location at x / D = 19.0 . 
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Fig. 4.14e Comparison between the predicted and measured turbulent fluctuation in the radial 

direction at x / D = 19.0 . 



Fig. 4. 14f Comparison between the predicted and measured turbulent fluctuation in azimuthal 

direction at measured location x / D = 1 9.0 . 
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Fig. 4. 15a Comparison between the predicted and measured mean axial velocity profiles at 

measured location x / D = 25.1 . 



Fig. 4.15b Comparison between the predicted and measured mean azimuthal velocity profiles at 

measured location x / D = 25.7. 
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Fig. 4.15c Comparison between the predicted and measured turbulent kinetic energy at measured 

location x! D- 25.1 . 



Fig. 4. 1 5d Comparison between the predicted and measured turbulent fluctuation in the axial 

direction at measured location x / D = 25.7. 
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Fig. 4.15e Comparison between the predicted and measured turbulent fluctuation in the radial 

direction at measured location x / D = 25.7. 



Fig. 4.15f Comparison between the predicted and measured turbulent fluctuation in the azimuthal 

direction at * / D = 25.7 . 
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Fig. 4. 16a Comparison between the predicted and measured mean axial velocity profiles at 

measured location x / D = 32.4. 



Fig. 4. 16b Comparison between the predicted and measured mean azimuthal velocity profiles at 

measured location x / D = 32.4. 
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Fig. 4.16c Comparison between the predicted and measured turbulent kinetic energy at measured 

location x / D = 32.4. 



Fig. 4. 16d Comparison between the predicted and measured turbulent fluctuation in the axial 

direction at measured location x / D = 32.4 . 


145 




Fig. 4.16e Comparison between the predicted and measured turbulent fluctuation in the radial 

direction at measured location x / D = 32.4 . 



Fig. 4. 16f Comparison between the predicted and measured turbulent fluctuation in azimuthal 

direction at measured location x / D = 32.4. 
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Fig. 4.17a Comparison between the predicted and measured mean axial velocity profiles at 

measured location x / D = 39.0. 



l-2r/D 


Fig. 4.17b Comparison between the predicted and measured mean azimuthal velocity profiles at 

measured location x/ D = 39.0. 
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Fig. 4.17c Comparison between the predicted and measured turbulent kinetic energy at measured 

location x / D = 39.0. 



Fig. 4.17d Comparison between the predicted and measured turbulent fluctuation in the axial 

direction at measured location x / D = 39.0. 
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Fig. 4.17e Comparison between the predicted and measured turbulent fluctuation in the radial 

direction at measured location x / D = 39.0. 



Fig. 4. 1 7f Comparison between the predicted and measured turbulent fluctuation in azimuthal 

direction at measured location x / D = 39.0. 
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Fig. 4. 19a Comparison between predicted and measured axial mean velocity at x / D = 2. 7 and 
28.5 ( U 0 is the axial mean velocity at the pipe centerline). 



Fig. 4.19b Comparison between predicted and measured azimuthal mean velocity at x / D = 2.7 
and 28.5 ( W 0 is the azimuthal mean velocity at the pipe wall). 
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Fig. 4.20a Comparison between predicted and measured turbulent fluctuation in the axial direction 

at x / D — 2.7. 



Fig. 4.20b Comparison between predicted and measured turbulent fluctuation in the axial direction 

at x / D = 28. 5 . 
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Fig. 4.21a Comparison between predicted and measured turbulent fluctuation in the radial direction 

at x / D = 2.7. 



Fig. 4.21b Comparison between predicted and measured turbulent fluctuation in the radial direction 

at x / D = 28. 5 . 
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Fig. 4.22a Comparison between predicted and measured turbulent fluctuation in the azimuthal 

direction at x / D = 2.7. 



Fig. 4.22b Comparison between predicted and measured turbulent fluctuation in the azimuthal 

direction at x / D — 28. 5 . 
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Fig. 4.23a Comparison between predicted and measured turbulent shear stress uv at x / D = 2.7. 



Fig. 4.23b Comparison between predicted and measured turbulent shear stress mv at x / D = 28.5. 


155 




Fig. 4.24a Comparison between predicted and measured turbulent shear stress vw at x / D = 2.7. 



Fig. 4.24b Comparison between predicted and measured turbulent shear stress vw at 

x I D = 28.5. 
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Fig. 4.25a Comparison between predicted and measured turbulent shear stress uw at jc / D = 2.7. 



Fig. 4.25b Comparison between predicted and measured turbulent shear stress uw at 

x/D = 28.5. 
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Fig. 4.26a Comparison between the predicted and measured axial mean velocity at four different 

axial locations. 



Fig. 4.26b Comparison between the predicted and measured azimuthal mean velocity at four 

different axial locations. 
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Fig. 4.26c Comparison between the predicted and measured turbulent kinetic energy at four 

different axial locations. 
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Fig. 4.27a Comparison between the predicted and measured rotation effects on axial mean velocity 

for different rotation numbers at x / D = 28.5 . 



Fig. 4.27b Comparison between the predicted and measured rotation effects on azimuthal mean 
velocity for different rotation numbers at x / D - 28. 5 . 
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Fig. 4.28 Comparison of the predicted rotation effects on turbulent kinetic energy at various axial 

locations. 
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Fig. 4.29 Sketch of coordinate system and computational domain for a fully developed flow in a 

rotating channel. 
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(a) 



Fig. 4.30 Comparison between predicted and LES mean velocity profiles in wall coordinates on (a) 
suction side (U + = U / u Ts , y + = (1 - y)u Ts / v ) and (b) pressure side ( t/ + =U / « r/; , 

y + = (1 +y)Urp / v) for /?o = 0.068. 
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(a) 



(b) 
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(C) 



(d) 


Fig. 4.3 1 Comparison between model predictions and LES data across the channel for (a) mean 
velocity; (b) turbulent shear stress; (c) turbulent kinetic energy; and (d) turbulent fluctuations. 
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Fig. 4.32 Comparison between predicted and DNS mean velocity in wall coordinates 

y + = (1 ± y)u r / 1> for Ro = 0. 
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(c) 



(d) 

Fig. 4.33 Comparison between model predictions and DNS data across the channel for (a) mean 
velocity; (b) turbulent shear stress; (c) turbulent kinetic energy; and (d) turbulent fluctuations. 
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(b) 

Fig. 4.34 Comparison between predicted and DNS mean velocity profiles in wall coordinates on 
(a) suction side (U + = U / u Ts , y + = (1 + y )u zs / v) and (b) pressure side 
(U + = U / u^, y + = (1 - y)Urp / v ) f° r R° = 0.01 . 
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(C) 



(d) 

Fig. 4.35 Comparison between model predictions and DNS data across the channel for (a) mean 
velocity; (b) turbulent shear stress; (c) turbulent kinetic energy; and (d) turbulent fluctuations. 
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(a) 



(b) 

Fig. 4.36 Comparison between predicted and DNS mean velocity profiles in wall coordinates on 

(a) suction side and (b) pressure side for Ro = 0.05 . 
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(d) 

Fig. 4.37 Comparison between model predictions and DNS data across the channel for (a) mean 
velocity; (b) turbulent shear stress; (c) turbulent kinetic energy; and (d) turbulent fluctuations. 
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(a) 



(b) 


Fig. 4.38 Comparison between predicted and DNS mean velocity profiles in wall coordinates on 
(a) suction side and (b) pressure side for Ro = 0. 10. 
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(a) 



(b) 
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(c) 



(d) 

Fig. 4.39 Comparison between model predictions and DNS data across the channel for (a) mean 
velocity; (b) turbulent shear stress; (c) turbulent kinetic energy; and (d) turbulent fluctuations 
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(a) 



(b) 

Fig. 4.40 Comparison between predicted and DNS mean velocity profiles in wall coordinates on 
(a) suction side and (b) pressure side for Ro = 0.15. 
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(d) 

Fig. 4.41 Comparison between model predictions and DNS data across the channel for (a) mean 
velocity; (b) turbulent shear stress; (c) turbulent kinetic energy; and (d) turbulent fluctuations. 
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(a) 



(b) 

Fig. 4.42 Comparison between predicted and DNS mean velocity profiles in wall coordinates on 
(a) suction side and (b) pressure side for Ro = 0.20. 
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(C) 



(d) 

Fig. 4.43 Comparison between model predictions and DNS data across the channel for (a) mean 
velocity; (b) turbulent shear stress; (c) turbulent kinetic energy; and (d) turbulent fluctuations. 
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(a) 



(b) 

Pig 4 44 Comparison between predicted and DNS mean velocity profiles in wall coordinates on 
(a) suction side and (b) pressure side for Ro = 0. 50 . 
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(C) 



(d) 

Fig. 4.45 Comparison between model predictions and DNS data across the channel for (a) mean 
velocity; (b) turbulent shear stress; (c) turbulent kinetic energy; and (d) turbulent fluctuations. 
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(a) 



(b) 

Fig. 4.46 Predicted mean velocity profiles for different rotation numbers in wall coordinates on (a) 

suction side and (b) pressure side. 
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Fig. 4.47 Comparison between predicted and DNS wall friction velocities u xp / u x and u TS / u T as 

function of Ro . 
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5 . VALIDATION AGAINST TURBULENT FLOWS WITH 
COMPLEX GEOMETRIES 


5.1 Introduction 

The NNWRS model is developed for complex turbulent flows, in particular for flows with 
complex geometry. Two-dimensional flow over a backward-facing step and three-dimensional 
flow in a square duct are typical complex flows resulting from the presence of multiple walls. 
Complex flow phenomena such as recirculation, reattachment, and boundary layer re-development 
to a fully-developed state are observed in the plane backward-facing step flow (Fig. 5.1). These 
flow features all result from the existence of the step. In a straight square duct (Fig. 5.10), the 
local flow structure is dominated by a transverse mean flow commonly known as secondary flow 
of the second kind induced by the presence of the corners. The mean transverse secondary flow 
consists of eight streamwise vortices, two counter-rotating in each comer, with the flow toward the 
corners from the duct center along the corner bisector, and toward the duct center along the 
bounding wall and wall bisector. In the following sections, two-dimensional backward-facing step 
flow and three-dimensional flow in a square duct are used as the testing cases to verify the 
NNWRS model's ability to predict flows with complex geometry. 

5.2 Two-Dimensional Flow Over a Backward-Facing Step 

5.2.1 Background Numerous attempts have been made to model two-dimensional 

backward-facing step flow by using various Reynolds stress models. Errors up to 50% - 100% 
between the model predictions of skin friction coefficient Cy = 2r w / ( pU %) and available data are 

often observed (So et al. 1988; Ko and Durbin 1994; Lien and Leschziner 1994). Here, U 0 is the 
bulk mean velocity at the entrance, p is fluid density and T w is the wall shear stress. One finding 
in the prediction of Cy by near-wall second-order models employing wall normals in the near-wall 
corrections is that the models tend to overshoot Cf in the adverse pressure gradient region after 

the reattachment point. However, no such overshooting is observed in predictions given by those 
high-Reynolds-number models invoking wall functions to satisfy the boundary conditions (So et 
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al. 1988) and Durbin's Reynolds stress model (1993). In Durbin's model, redistributive terms in 
the Reynolds stress equations are modeled by an elliptic relaxation equation to represent strongly 
non-homogeneous effects produced by the presence of walls, and no damping functions are 

needed (hereafter Durbin's model is referred to as the elliptic-relaxation model). Thus, the 
overshooting of Cy persisting in the near-wall Reynolds models is speculated to be caused by the 

use of wall normals in the near-wall correction terms. Since no wall normal is used in the 
NNWRS model, the comparison between the NNWRS model predictions of Cy and those given 

by the NWRS model will shed light on whether the use of wall normals is responsible for the 
overshooting of Cy in the near-wall turbulence modeling. 

5.2.2 Numerical implementations and boundary conditions The TEACH code and 
corresponding numerical procedures discussed in Section 4.2.2 are adopted for the current two- 
dimensional numerical computations in Cartesian coordinates. Only half of the expanded channel 
is considered since the flow is symmetric about the centerline. The computations are conducted in 
a domain 50hx6h shown in Fig. 5.2 with a grid distribution 91x81 in the streamwise and 
transverse directions, respectively. The grid is clustered at the inlet and the near-wall regions. 
Finer grid 131x151 was tested to give essentially identical results. Thus, grid 91x81 is used in 
the final computations. The convergence criterion is again given by (4-12). 

The model equations for two-dimensional backward-facing step flow are elliptic-type 
equations. To solve these equations in the (x, y)-plane, boundary conditions at the inlet and outlet 
of the computational domain, along the centerline of the expanded channel and the step walls are 
required. The mean velocity and the turbulent statistics at the inlet ( x = 0) are provided by the 
available DNS data. Negligible diffusion boundary conditions are specified at the outlet 
(x = 50 h), i.e. 

<9 2 (t/, V, uu, w, ww, mv, e) _ q , 5 _j, 

dx 2 

Along the step walls all the variables are zero according to the no-slip condition except the 
dissipation rate, which is given by 
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e = 2v\ 


OVP 2 

v dx J 


(5-2) 


along the vertical wall ( x = 0) and by 

OVP 2 


£ - 2v\ 


dy 


\ J 


(5-3) 


along the horizontal wall (y = 0). Symmetry conditions are applied at the centerline (y = 6 h) 
which yield 

d(U, MM, VV, WW, £) 


dy 


= 0 , and V = uv = 0 


(5-4) 


5.2.3 Cases studied and comparisons with data Both numerical and experimental study has 
been carried out for the same backward-facing step flow (Le et al. 1993; Jovic and Driver 1993) 
with Re h = U Q h / v = 5, 100 and the expansion ratio 8 / h = 1.2, where 8 is the boundary layer 
thickness and h is the step height. Detailed distributions of skin friction coefficient C f , wall 
pressure coefficient Cp = 2(P - P a ) / p(J^ (where P is the wall static pressure along the 
streamwise direction and P 0 is the wall static pressure at x = 0), mean velocity U , kinetic energy 
k, and turbulent stresses u rms , v rms , w rms , and uv at several downstream locations were given. 
The same case is used to test our models. 

Fig. 5.3 shows the flow pattern predicted by the NNWRS model. As a result of flow 
separation, immediately after the step, a primary cell and a corner cell rotating in opposite 
directions are generated. Further downstream, the flow develops into a plate boundary layer flow. 
The NWRS model predicts a similar flow pattern (no shown). Predicted friction coefficient C f , 
DNS results and measurements are compared in Fig. 5.4. The length between two zero crossing 
points of Cf is the recirculation region in the streamwise direction. The first zero crossing point 
signifies the beginning of the reverse flow. The comer cell locates between the step comer and the 
first zero crossing point of Cf. The NWRS model predicts the length of the corner cell is about 
2 h, which agrees with the DNS prediction, whereas the NNWRS model gives a slightly smaller 
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length. Durbin's elliptic-relaxation model predicts the smallest comer cell length among the three 
turbulent models. The second zero crossing point of C f is the mean reattachment point defined as 
the location of zero wall shear stress. The DNS and measurements show the flow has a 
reattachment length of 6 h, which is also predicted by both the NNWRS and NWRS models. The 
elliptic-relaxation model, however, overpredict the reattachment length, which is 6.6 h. The 
maximum friction coefficient (ICyl) in the reverse flow region is predicted well by the NNWRS 
and NWRS models, but it is underpredicted by the elliptic-relaxation model. After the reattachment 
point, the elliptic-relaxation model follows the DNS data and measurement quite well, although 
underpredicting Cy slightly. The overshooting of Cy in Reynolds-stress models mentioned by So 
et al. (1988) is also observed in the NNWRS model. This indicates that wall normals used in the 
near-wall corrections in the NWRS model are not the culprit for the overshooting. Further 
investigation is needed to find out the cause for this peculiar overshooting behavior in the Reynolds 
stress models. 

Fig. 5.5 compare the wall pressure coefficient C P . Model predictions and data all show the 
wall pressure along the horizontal step wall first slowly decreases from P a (indicating a favorable 
pressure gradient), and then rapid increases at some location in the recirculation region (indicate an 

adverse pressure gradient). It continues to increase and approaches to a constant as the boundary 
layer redevelops in the downstream. Contrary to their prediction of Cy, the predictions of C P 

given by the NNWRS and NWRS models agree very well with both DNS data and measurements, 
whereas the elliptic-relaxation model underpredicts C P in adverse pressure gradient region. 

Fig. 5.6 shows the comparison of normalized mean velocity U + = U / « T in the wall 
coordinate y + = u r y / l» at three downstream locations after the reattachment point: x/ h = 10, 15, 
and 19. The log-law and linear law are also included for comparison in the figure. The non- 
dimensional mean velocity profiles at these locations are predicted well by the elliptic-relaxation 
model, but underestimated by the NNWRS and NWRS models. This is consistent with the model 
predictions of the friction coefficient Cy after the reattachment (Fig. 5.4). The DNS data and 
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measurements show that the log region in the velocity profile at location x!h= 10 is very small 
and the intercept with U + - axis is also much smaller than 5.4 (Fig. 5.6a). As the flow develops 
further downstream, the length of log-law region increases and U + becomes closer to the log law 
(cf. Fig. 5.6a with Figs. 5.6b and 5.6c). The deviation of the velocity profile U + from the log- 
law may be due to the strong stream wise adverse pressure gradient (Nagano et al. 1991). It may 
also result from the non-equilibrium effects persistent after the reattachment. Further downstream, 
the adverse pressure gradient decreases, and the mean velocity profile U + becomes closer to the 
log-law. 

The performance of the models are further examined by comparing the model predictions of 
the streamwise mean velocity U and turbulent statistics (including kinetic energy k, normal 
stresses u rms , v rms , w rms > ar *d shear stress wv) with the DNS data and measurements in the global 
coordinate at locations x / h = 4, 6, and 10 (Figs. 5. 7-5. 9). Note that the location of the y axis in 
these figures is shifted for different measured locations. Fig. 5.7 plots the mean velocity profile 
U / U a versus y / h at x / h = 4, 6, and 10. It shows that U experiences rapid changes in the 
near-wall region and approaches a constant as y / h > 2. Unlike their predictions of U + j n the wall 
coordinate y + , the NNWRS and NWRS models give better predictions of U / U (} in the global 
coordinate than the elliptic-relaxation model does (cf. Fig. 5.7 with 5.6). In the near-wall region, 
all three models predict a negative U at x / h = 4, which agrees with the data. This confirms that 
location x / h = 4 indeed is inside the reverse flow region determined from the friction coefficient 
Cf (see Fig. 5.4). At x / h = 6, the elliptic-relaxation model predicts a small negative U in the 

inner layer, whereas the NNWRS and NWRS models predict a small positive U there. This 
indicates that the reverse flow region ends and reattachment begins near x! h- 6. After the 
reattachment point, at x / h = 10, all the models predict a positive U across the expanded channel; 
the NWRS and NNWRS models follow the data better than the elliptic-relaxation model does. 

The predicted turbulent kinetic energy k is compared with the DNS data in Fig. 5.8. Note 
that Jovic and Driver (1993) did not report the measured normal turbulent fluctuation w rms , so the 

measured kinetic energy k is not available. The mixing layer type structure in the flow separation 


193 


region behind the step significant enhances the turbulent production (Jovic and Driver 1993), 
resulting in high kinetic energy peak. As flow moves to the downstream, turbulence decays and 
the peak of the kinetic energy decreases. Compared with the DNS data, the NNWRS and NWRS 
models give pretty good predictions of k , although the NNWRS model slightly overpredicts the 
peak of k at x / h = 6 and 10, and the NWRS model underestimates the peak at x / h = 4. The 
elliptic-relaxation model underpredicts k at all three locations. The comparisons of turbulent 
normal stresses u rms , v rms , w rms and shear stress wv at the same locations are shown in Fig. 5.9. 

Overall, the elliptic-relaxation model given better predictions of turbulent stresses in the near wall 
region than the NNWRS and NWRS models, but the NNWRS and NWRS models outperform the 
elliptic-relaxation model in the outer region. The peak value of u rms is overestimated by the 
NNWRS and NWRS models, whereas the peak value of w ms is underpredicted by the elliptic- 
relaxation model. 

5.3 Three-Dimensional Flow in a Square Duct 

5.3.1 Background Turbulent flow in a square duct is characterized by the existence of 
secondary flow of the second kind (as classified by Prandtl) in the plane perpendicular to the 
streamwise direction. This kind of secondary flow in a non-circular straight duct is created by 
turbulent motion. Although relatively weak (2-3% of the streamwise bulk velocity), its effects on 
wall shear stress distribution, heat transfer rates, or transport of passive tracers are quite significant 
(Demuren 1990). Being of considerable engineering interests, turbulent flow in a straight non- 
circular duct has been the subject of many experimental and numerical investigations (Demuren 
and Rodi 1984). Systematic measurements of flow in a non-circular straight duct have been 
carried out by Gessner's group (Po 1975; Lund 1977; Eppich 1982; Gessner and Emery 1981). 
The farthest measured location in their experiments away from the entrance is 84 D h , where D h is 

the hydraulic diameter of the duct. According to the measurements, flow at this location is already 
fully developed. Similar experiments were also conducted in a duct with a length shorter than 
84 D h (Melling and Whitelaw 1976; Launder and Ying 1972). Although Brundett and Baines 
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(1964) carried out measurements as far as 260 D h from the entrance, where the flow is fully 

developed, the inaccurate method used to measure the turbulent stresses made their measurements 
less reliable. 


Direct numerical simulation provides an alternative approach to study the fully developed 
square duct flow at low Reynolds numbers. Recently, Huser and Biringen (1993) and Gavrilakis 
(1992) carried out direct numerical simulation for square duct flow at Reynolds number 
Re^ = U b D h / v = 10,320 {U b is the bulk velocity in the streamwise direction) and Re^ = 4,410. 

respectively. Although the Reynolds numbers in the DNS are much lower than those in the 
experiments mentioned above, it appears that the mechanisms driving the comer secondary flow at 
disparate Reynolds numbers are similar. Huser and Biringen (1993) showed that the low- 
Reynolds -number effects are manifested in the reduction of the distance from the vortex centers to 
the corner and the reduction of the secondary flow near the wall bisector compared to high- 
Reynolds-number experiments. Gavrilakis (1992) also found out that viscous effects are quite 
important in the transportation of mean vorticity in low-Reynolds-number flow, whereas they are 
important only in the region very close to the corner in high-Reynolds-number flow (Demuren and 
Rodi 1984). 


The origin of the secondary motion can be identified through analyzing the streamwise 
vorticity development along the duct. The mechanisms which cause streamwise vorticity to 
develop in turbulent flow along a corner are responsible for the initiation of the secondary flow. 
Following Prandtl's (1926) original idea on the secondary motion, Einstein and Li (1958) carried 
out a rigorous analysis to show that the gradients of Reynolds stresses are the actual source of the 
secondary motion. For a fully developed incompressible flow, the mean streamwise vorticity Q v 
is governed by 


V 


dQ^ 

dy 


+ W 


dz 



d 2 

l * 2 


d 2 \ 
By 1 ) 


vw+ v 




d 2 n x 

dy 2 



(5-5) 


where 
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(5-6) 


_dV_ dW 
x ~ dz dy 

In (5-5), the terms on the left-hand side represent the convection of the streamwise mean vorticity 
in the transverse plane; the first two terms on the right-hand side express the influence of the 
turbulent stresses on the production or destruction of the streamwise vorticity; and the last term on 
the right-hand side is the viscous-damping term. The two turbulent-stress terms are found to be 
the dominant terms, and it is the difference between these two relatively large terms that drives the 
weak secondary motion (Demuren and Rodi 1984). This implies that modeling the secondary 
motion in a duct requires an anisotropic turbulence model. 

Many attempts have been made to calculate the flow in a non-circular straight duct by using 
algebraic models or Reynolds stress models. The algebraic models can be obtained from the 
Reynolds stress models by neglecting the convection and diffusion terms. Most computations 
were carried out by using algebraic turbulence models (Launder and Ying 1973; Gessner and 
Emery 1981; Gosman and Rapley 1980; Naot and Rodi 1982; Demuren and Rodi 1984). 
Although the algebraic models can give reasonable overall predictions of secondary motion by 
tuning the model parameters, the models themselves do not correctly reflect the real physical 
processes (Demuren and Rodi 1984). Among the few applications of the Reynolds stress models 
to flow in a non-circular duct, wall functions were used to satisfy the boundary conditions (Naot et 
al. 1974; Reece 1976; Launder and Li 1994). The correct modeling of near-wall turbulence is 
crucial to the reproduction of the secondary flow in wall-bounded duct flow. Therefore, square 
duct flow provides a severe case to test the correctness of turbulent modeling of the wall effects. 
In addition, since the flow is bounded by multiple walls and has two inhomogeneous directions, it 

is also an ideal case to show the advantages of the wall-independent NNWRS model. 

5.5.2 Cases studied The high-Reynolds-number case (Re b = 250,000) investigated by 
Gessener’s group in the laboratory and the low-Reynolds-number case ( Rc h = 10,320) studied by 
Huser and Biringen (1993) using DNS are selected to test our turbulence models for square duct 
flow. Through the continuous effort of Gessener's group, fairly complete experimental data are 
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available for three-dimensional turbulent flow in a square duct at Re^ = 250,000. These data 

include mean velocities in streamwise and transverse directions, wall shear stress, and the turbulent 
stresses along both wall and corner bisector at streamwise locations x / D h = 8, 40 and 84. At 

location x / D h = 84, the flow is fully developed (Gessener and Emery 1981 ). These reported data 
can be found in Gessener et al. (1979, 1991, 1993) and Gessner and Emery (1981). Hereafter this 
high-Reynolds-number case is referred to as Gessner and Emery’s (1981) case. In Huser and 
Biringen's (1993) case, the flow is fully developed, and detailed DNS data are available. These 
two cases together will demonstrate the advantages of wall-indepent NNWRS model for flows 
with multiple walls and should complement the turbulence model validations over a wide range of 
Reynolds numbers for flow in developing region as well as in fully developed state. 

5.3.3 Numerical procedure and boundary conditions The coordinate system and 
computational domain for a quarter of the square duct are shown in Fig. 5.10. Three-dimensional 
parabolic marching scheme of Patankar and Spalding (1972) is adopted to solve the governing 
equations. This numerical scheme assumes that the streamwise mean velocity is determined by the 
averaged pressure in the plane perpendicular to the streamwise axis, and iterations are therefore 
carried out in the transverse plane. At fixed streamwise location ( jc ), the discretization procedure 
is similar to that described in Section 4 for swirling flow, and a line-by-line iterative scheme is 
used to solve the tridiagonal matrix in alternating directions. SIMPLEC algorithm (modified from 
SIMPLE) is adopted to link the pressure field with the velocity field. Fig. 5.1 1 shows the grid 
distribution of 91 x 91 in the lower left quadrant of the duct in the ( y, z)-plane. The grid points 
are clustered near the walls and 91 x91 is found to be fine enough for the cases studied. The 
marching step size increases progressively from 0.004 D h at the entrance (x = 0) to 0.04 D h at 
downstream location x / D h = 100. 

The inlet (x = 0) conditions determine the development of the flow in the square duct 
because the governing equations are of parabolic-type. At the inlet, the mean flow in the core 
region is assumed to be uniform with turbulence intensity I = 0.35%. Near the walls the flow 
field is obtained from the flat-plate boundary layer theory with relative boundary layer thickness 
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5 / a = 0.02 on the duct midplane ( y = a), where 8 is the boundary layer thickness and a is the 
length of the wall bisector (see Fig. 5.10). The secondary velocities V and W, and turbulent 
shear stress vw are set to zero at the inlet. All the variables are specified such that their iso- 
contours are parallel to the duct walls and have a square shape to give symmetric distributions 
about the corner bisector. No outlet conditions are needed for this parabolic problem. Symmetry 
conditions are used at the wall bisectors (y = a, z = a) , namely 


d(U, W, MM, VV, WW, MW, £) _ 

, and 

V — mv = vw = 0 , at y = a , 

(5-7a) 

dy 



d(U, V, MM, VV, WW, MV, £) _ 

, and 

W = mw = vw = 0 , at z = a , 

(5-7b) 

dz 




and no-slip boundary conditions require all the variables be zero at the walls (y = a, z = a) except 
for the dissipation rate, which is specified as 

e = 2v\ " I , at y = 0 , (5-8a) 


dy 


£ - 2v\ 




dz 


, at z = 0 


(5-8b) 


5. 3. 4 Results and discussion The model-data comparisons are presented first for the case of 
Huser and Biringen (1993) for fully developed flow at low Reynolds number and then for the case 
of Gessner and Emery (1981) for developing flow at high Reynolds number. Model predictions 
for the low Reynolds number case show that the flow at x / D h = 84 already reaches the fully 
developed state since no more change is observed further downstream. Therefore, the model 
results at jc / D h = 84 are used to compare with the DNS data given by Huser and Biringen (1993) 
for the fully developed flow. 

Fig. 5.12 displays the secondary flow velocity vectors predicted by the NNWRS model, 
revealing two streamwise, counter-rotating, and symmetric (with respect to corner bisector) 
vortices in the corner. The flow on the transverse plane moves towards the corner along the comer 
bisector ( y = z ), and by virtue of mass continuity, it moves away from the comer along the walls 
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( y = 0 and z = 0 ) and towards the center of the duct along the wall bisectors (y = a and z = a). 
This predicted flow pattern is very similar to that shown in Huser and Biringen (1993), with the 
centers of the vortices shifted slightly away from the corner. The vortex centers in Fig. 5.12 are 
located at (y ~ 0.25a, z ~ 0.55a) and (y = 0.55a, z ~ 0.25a), whereas their counterparts predicted 
by the DNS are located at ( y = 0. 2a, z = 0. 4a ) and ( y = 0.4a, z ~ 0. 2a ) 

Fig. 5.13 shows the isotachs (lines of constant streamwise mean velocity in ( y, z)-plane) of 
U / U c ( U c is the mean velocity at the centerline of the duct) predicted by the NNWRS model. 
The isotachs are bent toward the corner with a milder curvature near the corner bisector compared 
with those isotachs shown in Huser and Biringen (1993); the predicted isotachs close to the walls 
are also flatter. Similar secondary flow pattern predicted by the NWRS model is observed (not 
shown). 

The comparison between the predicted normalized wall shear stress t w / f w ( t h , is the 
average value of t w over the walls) and the DNS data along one wall is shown in Fig. 5. 14. The 
T w ! given by the NWRS model is much overpredicted near the corner, but agrees reasonably 
well with the DNS data away from the corner. On the other hand, the r w / r M , predicted by the 

NNWRS model follows the DNS data quite well from the corner to the location of the first peak in 
the DNS data; further away from the corner, the NNWRS model slightly overpredicts the wall 
shear stress. The predicted mean streamwise velocity profiles U along the wall bisector ( y = a) 

are plotted in Fig. 5. 15 and compared with the DNS data. Both models give very good predictions 
of the streamwise mean velocity U . 

In Fig. 5.16, the predicted secondary velocity V along the wall and corner bisector are 
compared with the DNS data. Both models overestimate V along the wall bisector except in the 
region near the center of the duct (Fig. 5.17a); the results given by the NNWRS model are closer 
to the data compared with those predicted by the NWRS model. Along the corner bisector (Fig. 
5.17b), although the NNWRS model follows the data near the corner and the center quite well, it 
predicts a secondary velocity profile with a much smaller peak. On the other hand, the NWRS 
model predicts a secondary velocity profile with a similar shape as shown by the DNS data, but the 
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peak is overpredicted by the NWRS model and its location is shifted toward the corner. The 
distributions of normal stress difference (w 2 - v 2 ) along the wall bisector are shown in Fig. 5.17. 
Compared with the DNS data, the NNWRS model gives a very good prediction of the normal 
stress difference, whereas the NWRS model much overestimates the difference near the walls. 

Fig. 5.18 shows the development of streamwise vortices in the lower left quadrant of the 
duct predicted by the NNWRS model for the high-Reynolds-number ( Re b = 250,000) case. The 
comparisons of secondary flow patterns at locations xl D ^ = 8, 40 and 84 indicate that the 
streamwise vortices in the developing flow are initially generated near the corner; as flow develops 
along the straight duct, they move away from the corner and spread out in the cross section until 
the flow becomes fully developed. At x/ D h = 84, the NNWRS model predicts the vortex centers 

are located at (y « 0.25a, z » 0.60a) and (y * 0.60a, z « 0.25a), which are further away from the 
comer compared with those predicted by the same model for the low-Reynolds-number flow (cf. 
Fig. 5.18c with Fig. 5.12). Fig. 5.19 plots the predicted and measured isotachs of the streamwise 
mean velocity U normalized by U c . The model predictions at x / D h = 8 (Fig. 5.19a), showing 

that the isotachs are bent toward the corner and are flat along the walls, are in very good agreement 
with the data. At jt / D h = 40 and 84, the predicted isotach contours are less distorted compared 

with the measured ones (Fig. 5.19b and 5.19c). 

The distributions of normalized wall shear stress x w / r w at three streamwise locations 

x/ D h = 8, 40, and 84 are shown in Fig. 5.20. Unlike the low-Reynolds-number case, the 
NWRS and NNWRS models give almost the same predictions of wall shear stress at all three 
locations for the high-Reynolds-number case (cf. Fig. 5.20 with Fig. 5.14). The predictions 
agree well with the measurements at x / D h = 84, although they are slightly overpredicted near the 

wall bisector. The rather flat predicted wall shear stress profiles result from the almost parallel 
distribution of isotach contours near the wall predicted by the models (Fig. 5.19). 

The variation of predicted centerline mean velocity in the streamwise direction is plotted in 
Fig. 5.21 and compared with the measurements. According to the experimental data, the mean 
velocity along the centerline first increases to the peak value, and then decreases and levels off to 
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an asymptotic value in the fully developed region. Both model results agree very well with the data 
in the developing region. In the shear interaction region near x / D h = 40, the NWRS model gives 
better predictions of U c than the NNWRS model does. The NNWRS model predicts a shorter 
developing region and a lower peak value of U c . As flow approaches to the fully developed state, 
the centerline velocity is also predicted well by the NWRS model, whereas underpredicted by the 
NNWRS model slightly. Fig. 5.22 shows the comparisons between the predicted and measured 
streamwise mean velocity U normalized by the bulk mean velocity U h along the wall and corner 
bisector at locations x / D h = 8, 40 and 84. The predictions given by both models are identical and 
in very good agreement with the data along both wall and comer bisector at all three locations. 

Comparisons between predicted and measured secondary flow velocity profiles along the 
wall and corner bisector at the same locations are shown in Fig. 5.23. Unlike the low-Reynolds- 
number case, both models give almost identical predictions of the secondary mean velocity for the 
high-Reynolds-number case (cf. Fig. 5.23 with Fig. 5.16). The secondary mean velocity is 
underpredicted by the models along the wall bisector at x / D h = 40 and 84. Along the corner 
bisector, the model-data agreement is good at x! D h - 40, but both models still underestimate V at 
x / D h = 84. The underestimation of V in the fully developed region is also observed when the 
secondary flow velocity is normalized by the local friction velocity u r , instead of the bulk velocity 
U b (Fig. 5.24). 

The comparisons between predicted and measured turbulent stress and kinetic energy 
distributions along the wall and corner bisector at three streamwise locations are shown in Figs. 
5.25-5.29. The NNWRS model predicts slightly higher turbulent stresses and kinetic energy than 
the NWRS model. Predicted turbulent stresses and kinetic energy along the wall bisector are in 
pretty good agreement with the experimental data. Along the corner bisector, the predictions of 
turbulent stresses and kinetic energy are generally overestimated. At jt / D h = 84, data and model 

comparisons are also made for the normal stress difference (w 2 - v 2 ) at locations z / a = 0. 1, 0.6 
and 1.0. Difference in the predictions given by the NWRS and NNWRS models is quite large in 
the inner layer, but diminishes away from the wall (y = 0). Compared with experimental data of 
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Eppich (1982), the normal stress difference is underestimated by the models in the near-wall 
region. 

5.4 Concluding Remarks 

In this section, the wall-independent NNWRS model has been applied to two turbulent flows 
with complex geometry: two-dimensional flow over a backward-facing step and three-dimensional 
flow in a straight square duct. The performance of the NNWRS model is shown to be as good as 
the wall-dependent NWRS model. The predictions of mean and turbulence flow fields in 
backward-facing step flow agree well with the DNS and experimental data. The NNWRS model is 
able to reproduce complex flow phenomena induced by the presence of the step, such as flow 
recirculation, reattachment, and boundary-layer redevelopment to fully developed state. The 
comparison of skin friction coefficient predicted by the two models indicates that the overshooting 
of the skin friction coefficient after the reattachment point prevalent in those near-wall Reynolds 
models using wall normals (including the NWRS model) is not caused by the wall normals 
adopted in the models. Further investigation is needed to find out the root cause for this 
overshooting behavior in Reynolds stress modeling. The NNWRS model is also shown to be able 
to capture the secondary motion induced by turbulence for low- and high-Reynolds-number flow 
in a square duct. It predicts that the vortex centers are closer to the corner in the low-Reynolds- 
number flow than in high-Reynolds number flow, which agrees with the finding of Huser and 
Biringen (1993). In high-Reynolds-number flow, both models give about the same predictions of 
the secondary flow velocity, whereas in low-Reynolds-number flow, the NWRS model predicts 
stronger secondary flow velocity along the wall and corner bisector compared with the NNWRS 
model. The larger difference between the model predictions for low-Reynolds-number flow seems 
to support the findings of Gavrilakis (1992) that viscous effects are more important in a low- 
Reynolds-number square duct flow. 
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Fig. 5. 1 Sketch of plane backward-facing step flow and coordinate system. 
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Fig. 5.4 Comparison between predicted and measured skin friction coefficient. 
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Fig. 5.5 Comparison between predicted and measured wall pressure coefficient. 
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Fig. 5.6 Comparison between predicted and measured mean velocity profiles U + in wall 
coordinate y + at different downstream locations after the the reattachment point: (a) x / h = 10, 

(b) x / h = 1 5 , and (c) x / h = 1 9 . 
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Fig. 5.7 Comparison between predicted and measured mean velocity profiles in global coordinate 

at different locations. 



Fig. 5.8 Comparison between predicted and measured turbulent kinetic energy at different 

locations. 
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Fig. 5.9 Comparison between predicted and measured turbulent stresses at different locations: (a) 

u rms , (b) v ms , (c) and (d) -uv . 
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Fig. 5.14 Comparison between predicted and DNS wall shear stress for fully developed square 

duct flow at Re = 10,320. 



a 


Fig. 5.15 Comparison between predicted and DNS axial mean velocity along wall bisector for for 

fully developed square duct flow at Re = 10,320. 
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(a) 



(b) 

Fig. 5.16 Comparison between predicted and DNS secondary flow velocity along: (a) wall 
bisector, and (b) comer bisector for fully developed square duct flow at Re = 10,320. 
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Fig. 5.17 Comparison between predicted and DNS normal stress difference along wall bisector for 
for fully developed square duct flow at Re = 10,320 . 


214 



215 






1 

0.8 

0.6 

a 

0.4 

0.2 

0 


Fig. 5. 18 Predicted secondary flow by the NNWRS model at: 
(a) x/D b = 8, (b) x/D h = 40 and (c) x/D h = 84. 
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Fig. 5.20 Comparison between predicted and measured wall shear stress for developing square 

duct flow at Re = 250,000. 



Fig. 5.21 Comparison between predicted and measured axial centerline velocity for developing 

square duct flow at Re = 250,000. 
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(a) 



(b) 

Fig. 5.22 Comparison between predicted and measured axial mean velocity along: (a) wall 
bisector, and (b) comer bisector for developing square duct flow at Re = 250,000. 
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(a) 



(b) 


Fig. 5.23 Comparison between predicted and measured secondary flow velocity along: (a) wall 
bisector, and (b) comer bisector for developing square duct flow at Re = 250,000. 
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Fig. 5.24 Comparison between predicted and measured secondary flow normalized by local 
friction velocity along wall bisector for developing square duct flow at Re = 250,000. 
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(a) 



(b) 

Fig. 5.25 Comparison between predicted and measured normal stress in axial direction along: (a) 
wall bisector, and (b) comer bisector for developing square duct flow at Re = 250,000. 
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Fig. 5.26 Comparison between predicted and measured normal stress in vertical direction along: 
(a) wall bisector, and (b) corner bisector for developing square duct flow at Re = 250,000. 
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Fig. 5.27 Comparison between predicted and measured normal stress in transverse direction along 
wall bisector for developing square duct flow at Re = 250,000. 
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(a) 



(b) 

Fig. 5.28 Comparison between predicted and measured turbulent kinetic energy along: (a) wall 
bisector, and (b) comer bisector for developing square duct flow at Re = 250,000. 
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(a) 



(b) 

Fig. 5.29 Comparison between predicted and measured turbulent shear stress along: (a) wall 
bisector, and (b) comer bisector for developing square duct flow at Re = 250,000. 
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Fig. 5.30 Comparison between predicted and measured normal stress difference at three axial 
locations for developing square duct flow at Re = 250,000. 
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6. CONCLUSIONS 


In this study, a new near-wall Reynolds stress model (NNWRS) without wall normals is 
developed based on the pressure-strain model of Speziale et al. (1991). With the aid of near-wall 
asymptotic analysis and the results of direct numerical simulation, wall-independent near-wall 
modifications are incorporated into the expression for the pressure-strain correlation, into the 
relation for the dissipation rate tensor, and into the modeled dissipation rate equation. A damping 
function is introduced to ensure that the near-wall modifications will not affect the flow field in the 
region far away from the walls. The asymptotically correct NNWRS model is the first near-wall 
Reynolds stress model without wall normals. For comparison, the formulation of the recent near- 
wall Reynolds stress model (NWRS) with wall normals developed by So et al. (1994a) is also 
presented. According to the asymptotic analysis, the NNWRS model gives more accurate 
predictions of those Reynolds stress components uncorrelated with the wall normal direction than 
the NWRS model, whereas the NWRS model gives more accurate predictions of those 
components correlated with the wall normal direction. With no wall-dependent variables used in 
the model and only one damping function, the NNWRS model is more general and flexible for 
turbulent wall-bounded flows than the wall-dependent NWRS model, which requires two damping 
functions for different Reynolds number range. 

The new model is applied to a wide variety of flows to verify its applicability. These flows 
range from relatively simple flows, such as fully developed channel/pipe flow, Couette flow, and 
boundary-layer flow with zero pressure gradient, to complex flows with swirl and rotation, such 
as swirling pipe flow, axially rotating pipe flow, and channel flow with spanwise rotation. The 
advantages of using the new proposed model for flows with complex geometry are demonstrated 
by flow over a backward-facing step and How in a square duct. The performance of the model in 
different type of flows is summarized as follows: 

(1) The NNWRS model predicts reasonably well the mean and turbulent flow fields of simple 

internal and external flows over a wide range of Reynolds numbers. The model is able to 


229 



capture the log-law region with a correct von Karman constant. Reynolds number effects on 
the mean flow and second-order statistics are reproduced well by the model. The main 
deficiency in the model is that it predicts less degree of turbulence anisotropy in the near-wall 
region compared with the NWRS model. 

(2) In flows with swirl and rotation, the flow complexity is caused by the streamline curvature, 
centrifugal force, and Coriolis force. The performance of the NNWRS model seems to 
depend on the relative importance of these factors. In swirling pipe flow, the NNWRS 
model gives reasonable predictions for flow in no reverse flow region, whereas the model 
predicts a much smaller reverse flow region than that observed in the experiment. In 
developing flow along a pipe rotating about its axis, the NNWRS model is able to capture the 
flow pattern subject to both destabilizing effect resulting from the mean circumferential shear 
strain and stabilizing effect due to the centrifugal force. In the fully developed channel flow 
with spanwise rotation, at low rotation rate, the model predictions agree with the data. But, 
as the rotation rate increases, the model underestimates the Coriolis effect. The NWRS 
model behaves similarly as the NNWRS model in these flows. 

(3) The advantages of using the NNWRS model for flows with complex geometry become more 
evident when it is applied to two-dimensional flow over a backward-facing step and three- 
dimensional flow in a straight square duct. The overall performance of the NNWRS model 
is found to be as good as the wall-dependent NWRS model. The model predictions of mean 
and turbulence flow fields in backward-facing step flow agree well with the experimental and 
DNS data. The NNWRS model is able to reproduce complex flow phenomena induced by 
the presence of the step, such as flow recirculation, reattachment, and boundary-layer 
redevelopment to a fully developed state. The comparison between the predictions given by 
the two models indicates that the overshooting of skin friction coefficient prevalent in near- 
wall Reynolds stress models is not caused by the use of wall normals. The NNWRS model 
is also able to capture the secondary motion induced by turbulence in low- and high- 
Reynolds-number square duct flow. The difference in the model predictions of low- and 
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high-Reynolds-number square duct flow seems to agree with the DNS findings by Gavrilakis 
(1992) and Huser and Biringen (1993). 

Overall, the predictions given by the NNWRS model agree reasonably well with available 
data from experiments, direct numerical simulation, or large eddy simulation for a wide range of 
Reynolds numbers, swirl numbers, and rotation numbers. Various complicated flow phenomena 
are essentially captured by the model. Moreover, the model performance is very close to that of the 
wall-dependent NWRS model. With further refinement, the wall-independent NNWRS model is 
expected to replace existing wall-dependent near-wall Reynolds stress models for wall-bounded 
turbulent flows soon. Modifications of the formulation of the near-wall corrections so that the 
NNWRS model can give more accurate predictions of those Reynolds stress components 
correlated with the wall normal direction are expected to improve the performance of the model 
significantly, which will be pursued in the near future. 
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APPENDIX A 


The Equations for the High-Reynolds-Number Models 

Speziale et al.'s (1991) high-Reynolds-number pressure-strain correlation formula was 
published in the following form 

n y = -( C 1 £ + C\P)b tJ + c 2 e{b lk b kj - In S^ + (C 3 - c 3 *n l/2 )^ 

2 N 


+ C 4 k\b ik S jk + b ik S jk - -b mn S mn Sij j + C 5 k(b ik W jk + b jk W ik ). 
The mean vorticity W tJ is defined as 


^ ij ®ij ^rnjr^-m 


(A-2) 


where a i)y is the anti-symmetric mean velocity gradient tensor, e mjl is the permutation tensor and 
Q m is the rotation rate vector of the non-inertial frame relative to an inertial frame. The definitions 
of the other terms in (A- 1 ) are 
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(A-3b) 

(A-3c) 

(A-3d) 
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and the coefficients are C, = 3.4, C 2 = 4.2, C 3 =0.8, C 4 = 1.25, C 5 =0.4, cf = 1.8, C 3 =1.3. 
Comparing (A- 1 ) with the LRR model, which can be written as. 
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(A-4) 


n f R = - C, eb tj - a, < P,j - | P5,J • ) - p, ( A, - 1 PS,J ) - 2 y, kS tJ , 


it can be noted that the first two terms in (A-l) are the return part and its nonlinear modification, 
respectively. The mean vorticity caused by rotation is separated from the anti-symmetric mean 
velocity gradient tensor and excluded from the transformation. Therefore, the terms need to be 
transformed are 


n ff G = (c 3 - c* 3 n ]/2 )ks lJ + c 4 k 

+ C 5 k(b ik O)j k + b jk co ik ) . 


b ik^jk fykSjk b mn S mn 5jj 


(A-5) 


The SSG model (A-l) was developed for plane homogenous turbulence by Speziale et al. 
(1991). Under this assumption the symmetric and anti-symmetric tensors S {j and 6)^ can be 
written as 


(A-6) 



r 0 

5 

0^ 


^0 

Q) 

0^ 

■% = 

s 

0 

0 

; = 

-O) 

0 

0 


1 ° 

0 

0y 


v0 

0 

0; 


and the anisotropy tensor is of the form 
0 A 


b tj = 


(bji b l2 

b }2 b 22 
0 0 b 


0 


33 ) 


(A-7) 


Accordingly, 


nf c = (C 3 -C 3 *n 1/2 )jt 


( 0 5 0^1 
5 0 0 

v 0 0 0 j 


+ C^k 


~b J2 S {bu + b 12 )S 0 

(b n +b 22 )S —b/ 2 S 0 


0 


0 


+ Czk 


2b I2 (0 {p 22 ~ b\ ] )co 0 ^ 

[p22 ~ b\ i )o) -2b I2 o} 0 

0 0 0 


;b/2^ 
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or 


n 


SSG 


' c II c, 2 0 'l 

Q2 ^22 0 


V 


0 0 c 


(A-8) 


33 J 


where 

Q 1 


jC 4 +2C 5 jb l2 kS , 


Q2 ~ [(Q _ C 3 n ' /2 ) + C 4 (b n + b 22 ) + C 5 (b 2 2 -b n )^kS , 


C 22 = | -Q-2C 5 \b I2 kS , 


C 33 - -—C 4 b 12 kS , 


(A-9a) 

(A-9b) 

(A-9c) 

(A-9d) 


To determine the coefficients in the final form (A-9), consider a thin shear layer flow where 

dU 


S= co = 


2 dy 


(A- 10 ) 


On substitution of (A- 10) and the anisotropy tensor definition (A-3a), (A-9) becomes 


C,, = - 
11 2 

( 

C] 2 = 


l r V dU 
^C 4 + C s j„v— , 


Q _ Q. + Cj ~ Qn 
12 4 4 


1/2 A 


2 


C 22 - 


f 

V 

If 
2 k. 


Q ! Q ~ £3 n 

6 4 


*nl/2\ 


W 


— du 

u + 

dy 

d U 
dy ’ 


r Q , c s | c, - c,’n ' <2 ~W at/ 

12 4 4 I ,9y 


C 4 - C 5 |WV 


<9E7 

<9y 


~ I r —dU 
Ca = --Q„v_ 


(A-l la) 


(A- lib) 


(A-l lc) 


(A-l Id) 


Similarly, last three terms in (A-4) can be written in the form 
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or 

'c], Cj2 0 ' 
nf c = c n c !2 o 
0 0 c\ 3 

\ 5i J 

where 

, 2 — dU 

C,, = (2a, -/?,)- mv— 



h 

2 






dU 

dy 


• 2 , n x — dU 

C22 ~ “ a ') uv ~fy 


■ 2 / n^-du 

C 33 - — j( a l + A) WV ^ 


By comparing (A- 11) with (A- 13), the following equations are obtained 

|( 2 a,-A) = tQc 4 + C 5 ' , 

9 \ f \ \ 

±( 2 p l -a l ) = -[-C 4 -C 5 ] , 


(A-12) 

(A- 13a) 
(A- 13b) 
(A- 13c) 
(A-13d) 

(A- 14a) 
(A- 14b) 
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(A- 14c) 


-(a, + P\) = -C 4 


a 1 L = ^4+£5 + i( c C 3*n 1/2 ) , 

*2 12 4 4 V 3 3 > 

^-f=f-f4( c 3- c 3n ,/2 ) . 
-T = -% + l( c ^ c ^ n]/2 ) • 


6 4’ 

Only three of the above equations are independent and undetermined coefficients, a } , /?, 
are solved from them, which are 


a, = 


.OiQ . a = Q^Q . *-%-±(c,-t$n '«). 


4 ' ' 1 4 ' ' ‘ 3 2 ' 

Substituting C 3 = 0.8, C 4 = 1.25, C 5 = 0.4, C 3 * = 1.3 into (A- 15), we have 


a, = 0.4125 ; )3, =0.2125 and y, = 0.01667 + . 


Therefore, the transformed SSG model is 

nf c = -(C,e+ c; p )b tj + - 5 n5 0 ) + c s a m (b ikemkj + V™*, )* 

with ctj =0.4125; /?, =0.2125 and y, =0.01667 and C, =3.4, C 2 =4.2, C 5 =0.4, 


(A-14d) 

(A-14e) 

(A- 140 
, and y, , 

(A- 15) 
(A- 16) 

(A- 17) 
:,* = 1.8. 
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APPENDIX B 


SUMMARY OF MODEL EQUATIONS 


B.l NNWRS Model 


The Reynolds-stress transport equation for the NNWRS model can be written as 


du t Uj _ d 

U i 


f duM, 


dx L dxi 


V- 


i-j 


V dx k ) 


du, 


du t 


+ D'J + [~ Ui u k - UjU k — 


dx k J " K dx k 


2£l k (UjU m e !kin + ) + Yljj £,y . 


where the modeled terms are 


Dl= 3 


lJ dxi. 


C r — 


du :U k 

u , u l — 
ox i 


du k u i 


dujUj ^ 


+ UjUl ~d^ + UkU ‘ dx 


1 ) 


Uij - -(C,£ + C, Pjby + c 2 e(b ik b kj - - n Sij ) + C 5 n m (b ik e mkj + b jk e mkl )lc 
nj = (C, £ + C, P)b,j - C 2 e(b lt b kj - ±n«<,)+ a\p v - \ PS kJ ) + Iffy , 


£ ij = - S ij £ ( 1 - fw.l ) + fw. 1 e 7j + £ ij - 


vv & 

£,j =-U iUj 


* 1 

£ ii = - 
J 2 


P d UjU . M.M | 


dXr. dXr 


U;Ui d dj t 

(V— ) 


k dx, dx b 


The dissipation rate equation is given by 


U, 


de 


dxi dx,. 


de 

ydx kJ 


de 


+ *“ (Ce 7“'"' lk, )+Ce '~k F ~ Cclfc T 


+ C f3 v — 

£i k 


d^fk 


\ 2 


V dx k J 


(B-l) 

(B-2) 

(B-3a) 

(B-3b) 

(B-4a) 

(B-4b) 

(B-4c) 

(B-5) 
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B.2 NWRS Model 


The Reynolds-stress transport equation for the NWRS model can be written as 


U, 


dujUj _ d 


dxu dxi 


duM 




\ V dx k ) 


dU ; 


dU: 


+ fyj + [~ u i u k ~ “j u k 


2.^1ki u jk t m^ikm jkm^ ^ ij ^ ij ■ 


(B-6) 


The modeled terms are given by 


D t = — 
lJ dx. 


C- 


dujU k 


V 


dx 


du k u. du^: 

“3 - + u k u l-^- ~ 
uX[ OX 1 


(B-7) 


n y = -(C]£ + Q P )b [J - + C 2 £(b ik b k j - -n5y ) + C 5 Q. m (b ik e mk j + bj k e mki )k 


- a, (Pq - | P<5„ ) - A ( D„ - 1 PS,J ) - 2( y, + n 1 1,2 )kS :j 

+ /», .n,“' + Ilf • 


(B-8a) 


where 


m = — 
y 3 


<9x ; 


i>- 


(9 m , u k 
dx 


1 ) 


+ (9x, 


c 

3 


^ (9m, M ^ 


V- 


j k 


dx 


1 ) 


£ ij =^ijS(l-f wA ) + f^ 


w £ UjUj + U i u k n k n j + UjU^n, + n,n } 
£,J ” k 


1 + 3u k uin k rii / 2k 
The dissipation rate equation is given by 


U k 


de 


dx k dx k 


de 

v l~ 

\ dx kJ 


dXj 


k - 

— i 
£ 


de , 

dx ; 



(B-8b) 

1 d ( duM , ^ 

3*.r *. r w * 

(B-8c) 


(B-9a) 


(B-9b) 

p-c £2 f + C , 

(B-lOa) 


where 
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£ = /, 


w,2 


( -2 ~ \ 

£ ~ £ ££ ' 

-L-P + M—-N — 
y k k k j 


£ = £ - 2 vk / v 2 , 

£ = £-2 v[d^fic / dyf , 

Other tensors and vectors appearing in the NWRS and NNWRS models are given by 


P U=~ 


du, 


J dUi 

u i u k n + u j u k -T— 


<?jr 




^ = 2P„ , 


blJ = 2k {U,Uj ~3 k6lj) ' 


r-r l L 


y 2 cfoy 


a,= 


<9C/, 

U.Ui, — — + U ,Ul 

dx, J k dx, 


n = (0,1,0) , 


e y* 


0, // any two of i, j, k are the same 
< 1, ifijk is an even permutation of 123 
-I, ifijk is an odd permutation of 123 


(B-lOb) 

(B- 10c) 
(B-lOd) 


(B- 1 1 a) 

(B-l lb) 
(B-llc) 
(b-l Id) 
(B-l le) 

(B-l 10 
(B-l lh) 

(B-l li) 
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B.3 Model Constants and Damping Functions 


The model constants and damping functions are listed in the table below. 


Constants/Functions 

NNWRS 

NWRS 

c, 

3.4 

3.4 

c 2 

4.2 

4.2 

c; 

1.8 

1.8 

C* 

U 3 

1.3 

1.3 

a, 

0.4125 

0.4125 

A 

0.2125 

0.2125 

7i 

0.01667 

0.01667 

a 

- 0.32 

- 0.29 

y 

0.072 

0.065 

c s 

0.11 

0.11 

C', 

1.5 

1.5 | 

C £2 

1.9 

1.83 

C e3 

2.95 


C 5 

0.4 

0.4 

C e 

0.12 

0.12 

L 


2.25 

M 


0.5 

N 


0.57 

fe 



f wl 

e -(Re,//50) 2 

ARe t /60)* ^-(Re,/200) 2 

fw2 


e -{R',/40) ! 
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where 


Re, 
A = 1 


V£ 

-(9/2)(b iJ b lj -2b IJ b ji b ki ) 


(B- 1 2) 
(B-13) 
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APPENDIX C 


SUMMARY OF GOVERNING EQUATIONS 


C.l Fully-Developed Channel Flows 

C.l A Mean flow equation 
d 2 U + duv + 1 


0 -- 


-+■ 


dy + dy* Re T 


(C-l) 


CA.2 NNWRS modeled equations 


0 = 


dy + 


( 


k* -t + ^ 




/ + C S — v 

e J 


k + — + duv 


du 2 . 

— - + 2 Q — mv — 
dy e dy 


- 2 “ v ^7 


-c, £*(/- + {|(«, - /„,«•) - §/>, + c;(/-/»x} 


- + <9tT 

MV o>y + 


+C 2 e + (1 - /„, )[b u 2 +b ]2 2 -^ny^{\- f w] )e + - /„, pr « 2 + 


1 

Qj 

( du 2 + ^ 

u 2 d 

(3k*) 



k + dy + 

& 

+ 

1 


(C-2) 


0 = 


<?y H 


V 


1 + '5C s —rV 2 
e* 


+ W + 1 


dy' 


C x e + {\-f wX )b 21 


-{§(«/ - /*/«*) - ~ 3 P> ' c ^( 7 - /./)^}«v + |p- 

+C 2 e + (1 - /„, )^ 22 2 + ^, 2 2 - |n] - 1(1 - )e + ~ f wl pr v 2+ 


1 

d 

' dv 2+ ^ 

v 2 d 

f dk + ) 


2 

dy + 

l* + J 

k + dy + 

Wj 



(C-3) 
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0 =- 


d V + 


fc + — ^o>w 2 
/ + C c — v 




By* 


-C,£ + (l-f wl )b u 


-{I (“r/y) + 1 A - c /*( 7 - /»/K} mv+ £ 


+C 2 e + (l-/ vvl )(&33 2 -^n]-|(l-/ vvl )e + -/.,^vv 2+ 


^w 2+A 




<?y 


VV 


2 " 3 ^ 


k + dy + 


\ J 


dy 


\”.r J 


0 =- 


dy' 


, ,~ r k + 2 + > l^«v 
1 + 2Cc— -v 
* £ + 


J 


, „ k + - + dv z 
— p + Q — mv — — 
dy e dy 


~2+dU* 


- V 


By* 


-c,£\l- f wl )b, 2 +(a,- f w ,a) 

-+ dU + ( 


+c,(l- -f wl )b l2 uv 


dy* 


•\~ + dU + n ~ + dU* 

J7 +P ‘ u ' 

r-f.'f+^n"* | k 


■ c] nm v. sv' 


By' 


^»+ 

+ C 2 £ + (l -fw\){b] 1 ^12 + ^12 ^ 22 ) _ fwl ~j^+ UV 


B f 


^ duv + ^ 

v^y 


uv 


f dk + ' 


k + dy + 


K*J 


(( 


0 = 


By* 


\\ 


\ + C c ^ 


de H 




- £ + - + dU + £ + , + 

~C F ]—UV - - -C £l —rr 

£] k + dy + k + 


+ C. 


c3 ' 


d4k+ 

v y 


(C-4) 


(C-5) 


(C-6) 
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C.1.3 NWRS modeled equations 


0 '= 



c + - + W 




+ 2C, 


-uv 


duv + 

w 


- 2uv 


dU + 

dy + 


+c,£*(/-/»,)(v+V-jn) 




vW 


A + + 5v 2 /2 


e+ -7( 7 -/^K 


(C-7) 


f? 


0= 


<?/l 


- + 3C 


k + -7 




<9v 2 




-{f (“, - /.,«') - - c;(i- f.,)b J2 }«v* 

+C i £*(/-/.,)(* 2; 2 +6„ i -j;7j 


A + +5v 2 /2 i 


(C-8) 


0=- 




1 + Cr—V 2 ' 




dw 


r+A 


-c, «*(/-/.>„ 




w 


-/wi 7 "-+■ g + --(/-/wK 

A: + 3v 2 /2 •* 


(C-9) 
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0 =- 


df 


f 2 k + —+^\duv + „ k + — + dv 2 

— + 2C, — v 


-> + A 


<?>> 


+ C c — MV 


+ 1 '-5 


dy + 


— v 


dtr 

dy + 


~C, £ + {l-f wl )b l2 + {a, - L'Cc'y'^r + P'U 2 * dU 


d/ 


df 


+ C]{l-f w ,)b I2 uv 


dW 

dy + 


Yi ~ fwi y + 2 77 


. c; „, n \+Bu+ 


dS 


+ C 2 £*(l- f wl ){b n b, 2 + b l2 b 22 ) - 2/, 


uv 


wi — -+ ' 

k + 3v /2 


(( 


0 -- 


dy + 

+ f»2 


\\ 


, + c, k -7 

e J 




dy + 


£ + — + dU + „ £ + £ + 


” c " r “ v W ~ c “ ** 


r T £ + - + dU + A/l £ 2 .,£ + £ +N 

L—uv — — + M N — - 


V 


r dy 


J 


(C-10) 


(C-ll) 


C. 1.4 Boundary conditions 
At y + =0 (wall) 


U + = 0, k + = 0, u 2 =0,v 2 =0, w 2 = 0, uv + = 0 


£ + =2 


f d~Jk 


V 


v dy* j 

At y + = Re r (centerline) 


dU + dk + du 2 dv 2 dw 2 — + 

17"°’ J7" 0 ’ W= 0 ’ W= 0 ' ' " v 


d£* 

dy* 


-0 


(C- 12) 


(C-13) 
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C.2 Fully-Developed Pipe Flow 


C.2.1 Mean flow equation 


0 = 


1 d 

r* dr* 


r* 


dU *' | 
dr + ) 




(C-14) 


C.2. 2 NNWRS modeled equations 

f f 

- r 1 + CV — v 

t, l £ J 


0= 


/ <? 


- + dr 


du 2 * + k * — + duv + ^ 

+ 2r Cc — mv 

5 e dr + 


dr* 


- 2uv 


J 


-C, £*(/- + {|(a, - /„,«•)■ -|j8, + c;(/- 


+ 

dr* 

dU* 


dr* 


1^2,. g. \ _i_ . e + ~2 + 


+ C 2e + (l-/ H , 1 )[fe 1 , 2 +^2 2 -3nJ--(l-/ vvl )e + -/ tvl ^ r « 



f 

— ~ + 'N 


1 d 

r + 

du 2 

w 2 1 <9 

r + dr + 

r 

dr + 

* + r + <?r + 


l 

J 



^ + <& + ^ 

v r 


0= 


>- + <?r + 


' ' uv ‘‘TlJ? 1 

1 + 3Cc — v 
4 e + 


v 


v 


y 


dr 4 


2 f-j+ 72 + ^ 2 2v" 


V - vv 


J 


3 ^ 


2 „ k*~ r dw 2 4 „ k 

r* e* dr 


4 _ k* — +/— + — + \ 

~ c v ( v ) 




— ♦ (?£/ + 
“ V <9r + 


+C 2 e + (1 - /„, )^ 22 2 + b ]2 2 - 1 n ) - 1(1 - /„, )e + - f w] ^ v 2+ 


1 d 

( -,+\ 
+ dv 2 

*- 

v 2+ 1 d 

f + dk + ) 

1 

+ 

-F 

dr + 

\ J 

k + r + dr + 

{ dr + J 


(C- 15) 


(C-16) 
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0 = 


1 d 


r* dr * 


l+C s —v 

\ £ J 


dr* 


2 /— + — +\ 2 v 2 

r +2 \ ) 3 r +2 


i 2 r k+ ~2*dw 2 , 4 r k+ ~*(~* — + \ 2 <9 (_ /fc + — ♦/— + — ^ 

V c v 5 r * + t tC *f w ( v _w rv" r 


■C, eV-f.,K - {|(«, |ft - c;(/- ^ 


+C 2 e" (1 - )( V n)-|(i - /„, )£* - /„ C„, 2+ 



f 

-T+ \ 

1 <9 

r + 

dw 2 

r + <9r + 

r 

dr + 


V 

J 


w 2 1 d ( + dk + ^ 


k + r + dr + 


dr 4 


v y 


0 = 


1 <9 


- + <9r + 


f £ + -7+^ ■* + 

1 + 2C S — v 2 


) 


duv + + k + — + dv 2 ^ 
— r + r c s —^uv — 

J 


d^ 


2 MV 

3 ^ 


l„k*— + dw 2 2 k* — + — + -+<9t/ + 

-Q— MV — — --^Q— MV w -V — T -C,e {l-f wl )b l2 


r* > £* dr* r* 2 ' £ 


dr 4 


+( a i-fwi*x*y ^r+^/W 2 ^r + C( 7 -/../K, MV+ ^r 


A n #y+ 

7 , - /„! 7* + y n in y — + C 2 e*{l- f wl ){b„ b, 2 + b l2 b 22 ) 


£ + — + 1 
/ji ; 

1 

^ + duv + ^ 

MV + 1 <9 

Ml/ 

£ + 2 

r + dr + 

**J 

& + r + <9r + 


( + dk + ^ 


V 


0 = 


1 d 


- + < 9 r + 


+ c e3 1 * 


r, 

<9e + ] 

1 4- C F — — V 

p + 

v £ ; 



V 


_ e*-+BU* „ £* s * 

-C £l — MV — ; Cl — T P 


fc + "' ( 9 r + £2 * + 


(C-17) 


(C- 18) 


(C- 19) 
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C.2.3 N'A'RS modeled equations 


0 = 


7 d 


' + dr + 


( ( 


k + ~Adu 2 + , _ + _ k + — + duv 




/ + C 5 -v 


ir' + 2rC s -uv Sr , 


— 2uv 


+ dU + 

dr + 


-C, e + {l- f wl )bj, + 1 ~ 3 (a, - f wl a ) - j/3, + C]{1- f wl ]b,^uv + ^ 


+C 2 £ + (7-/ M „)[ b „ 2 + b n 2 -^n 


-f wI u 2+ L— --{i-f w y 
k + 3v 2 / 2 3 


0 -- 


1 d 


r + dr + 


(r 




2 , r /c + --yv 2+ ^ 

- + 3C S — v —— 

3 V J< 9 r + y 


2 /— + — +\ 2 2v 2 

+2 (v w ) + 3 r +2 


2 „ k* - r+ dw 2 4 _ k 


r Cc — v 


r + £ + < 9 r + r + " J £ 


4 _ r— +/— + — +\ 

--TtQ-h. (v -W ) 


*\ 4 


-C, e + (l~ f wl )b 22 - -la, - /.,«)- -B, - C]{1 - f wl )b 22 


— + dW 


uv 


dr 4 


+C 2 £ + (7-/ K , / )( 7j 22 2 + ^/2 2 


~ 4 fwi y2 — .. ,(/ /»/) £ 


r+3v 2 /2 3 


0 = 


7 <9 


' + < 9 r + 


7 + Q~7 y2+ 




V 


< 9 r + 


+ 


2 ( 2 + 2 + \ ^ V 2 

r V -W 2 

r +2 1 / 3 r + 


2 ^ k + 2 + <9w 2 


r £ 


4 _ k* —j+(- 2 + ~ + \ 2d [ _ k + —+(—* — j + 

+ 7 C V W ( v " w ) + 797 ( C 'F" ( v " w 


< 9 r + 


-C, e + (7 - - { j(«, - /„,«*) + f ft " ^*(7 - /wKW 


(C-20) 


(C-21) 
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+c 2£ + {i-f w ,)[b }3 2 -- 3 n 


-fwl* 


k + + 3v 2 /2 3 




0 = 


1 d 


dr+ 


+ 2 , _ k + -r 


^ d uv* , + _ k + — + dv 


r\ — + 2C c — 


3 


£ ) 


, + r C, — -uv 
dr 


_2 +> \ ~ + 

2 uv 


e + dr + 


—Lr—w +dw2+ 2 r k+ ~ + ~ + ~ +dlJ 

r«v — — ; — uv w ~ 


r J e + dr 


+ +2 '-'5 

r e 




+(«, - /.,«v ‘ + p,u ! ' jp- + c,-( l - /„, )»„«,■ 


— + dLT_ 

dr + 


+^-n" 2 


dU 


k+ ~ df +Cl £+ ^~ */2 + */2 *22 ) 


(C-22) 


o /• + £ 

~2f w ,uv 


k + + 3 v 2 /2 


0 = 


1 d 


f 


r + dr 

+ fw2 




V c 


<?£ 




y<9r 4 y 


^ £ + — + <?£/ + e + e + 

-C ..—uv C, 

k + dr + E2 k + 


f r e + — + dU + 

L—uv 

k + dr + 


e 2 e + e + ^ 

+ M- N-— 

k + k + j 


(C-23) 


(C-24) 


C.2.4 Boundary conditions 
At r + = Re T (wall) 

U + =0, k + =0, u 2 =0, v 2+ = 0, w 2+ = 0, uv + - 0, 



d-^fk 

~) 


y 

J 


At r + = 0 (centerline) 


(C-25) 
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=0, 1^=0, d -fr=0. i Tv‘ = 0, 

<?r + <9r <9r + <9r <?r 


<?e + 

<?r + 


=0 


(C-26) 
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C.3 Fully-Developed Couette Flow 


C. 3. 1 Mean flow equation 

<? 2 (/ + dw + 

~ dy +2 " 

C.3. 2 NNWRS modeled equations 

The equations are the same as those given for the fully-developed channel flow. 


C.3. 3 NWRS modeled equations 

The equations are the same as those given for the fully-developed channel flow. 


C.3. 4 Boundary conditions 


At y + =0 (fixed wall) 


U + =0, k + =0, u 2+ =0, v” H 

- /MY 
“ U/J ■ 

At y + = 2Re T (moving wall) 


0, w 2 = 0 , uv = 0 


u + = u: pecified , k + =0, u 2+ - 0, v 2+ = 0, w 2+ = 0, uv + = 0 


£ + = A 


djk \ 2 

\Jy + ) 


(C-27) 


(C-28) 


(C-29) 
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C.4 Incompressible Boundary-Layer Flow 

C.4. 1 Mean flow equations 
dU dV „ 

<9x <9y <9y <9y [ dy 



C.4. 2 NNWRS modeled equations 


u^L +v ^± 

dx dy dy 




V 


v + Cc-V 


du 2 k — duv 

_- + 2Cj -„v— 


— 2 MV— — 
o>y 


+C 2 £(1 - /„, )(*i i 2 + i|2 2 - 1 n j - 1(1 - /„, )£ - t-u 2 


1 

d 

d? n 

i > 

M 2 <9 

r ^ 

-j \ 

2 

dy 


& <?y 

l *J 


u<hi +v i?L = ± 

dx dy dy 




v + 3 C s -v 2 
6 e 


— \ D ..2 


<9v 2 

<?y 


- C] £(l ~ f w \)b 22 


-{ 3 («, /wl 01 ) "jA — Q 0 _ /w] )^22 j MV 


— du 

dy 



+C 2 e(l-f wl )\b 22 2 +b n 2 --n -±{i- fwl ) e -f* v 


£~2 




dy 


v- 


dv 


dy 


dk 


\ 

cm 

k dy^ dyj 


U 


dw 2 

dx 


+ V 


dw 2 

dy 



_ k ~ 2 \ dw 
v + Cc — v 1 
e 


2 ^ 


J <?y 


- C]£(l ~ f w \)b 33 


(C-30a) 

(C-30b) 


(C-31) 


(C-32) 
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(C-33) 


_ {j(“l _ /v«l a ) + |A -C|‘(l-/ rf )*3 3 J„v|y 
+c 2 e(I - )(i, 3 2 - In] - 1(1 - /„, )e - 

_1 d v dw 2 W 2 <? r 

2 <?y J k dy\ dy J 


duv + duv _ d 
dx dy dy 


v + 2CV — 1 / 


A: 2 ^ 


T- + C s -«v~- 

dy e dy 


y2 dy C ' € ^ ^>>i2+(ai-/vvi« )v 2 -|y + AM 2 -|y 

,^* n / n*. ~ dU ( S * c 3 rrl/2^ du 

+C\{\~ f w \)b n uv— Y\ ~fw\Y +-r-n 1/2 k—~ 

dy y 2 J dy 


+ c 2 £ { l -fw]){bnb ]2 +b l2 b 22 )- f w ]juv 

k 

1 d f uv d ( dk ^ j 

— — U— — v — 

2 L ^ V dy J k dy{ dy )_ 


J T C/0 T 7 UO V „ 

U— + V—=— v + c, 


k-j^de 


C.4.3 NWRS modeled equations 


TT du 2 ,,du 2 d 

U — 1- V— — = — 

dx dy dy 


+ 2 C s -uv^- 
dy e dy 


duv | „ — dU 


-C\ ~ fw\ )b\ 1 + jj(«l ~fwl a )"A +C ’i*( 1_ /n'l)^ll| wv ^ 

+c 2 e(i - )(i>, 2 i + *i 2 2 - j n> 


(C-34) 


(C-35) 
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(C-36) 


fw\ ~2 £ 3O fw\) e > 

k + 3v 2 / 2 3 


^ + v 2£4 

(7 X <7 v <7 y 


0 \ 


2 k ~2\dv 2 

— D + 3CV — v 


e J <9;y 


-{|(«, - /wl« ) - f A - C 1* 0 - /wl )*ta}«V 


+C 2 e(l-/ H ,, 1 )(fo 2 2 2 +^i 2 2-^n) 


"^/wl" 0 £ 3^ /H>l) e * 

k + 3v l /2 3 


C,e(l-/»i)*22 

i-au 






v u+Cs 7 v J~77 


<9 jc <9y dy 

{| («1 - /wl a * ) + \\ A - cf ( (1 - . /wl ) )^3 3 ^ 


-Ci(l“/wl)^33 

}-du 


Jl 1 


+ C 2 £(1 _ /wl)(^33 


-/, 


w 


wl 


k + 3v 2 /2 3 


£ ~ (l /wl)^ ’ 


ri <? MV „^ uv 9 

l, -57 +v 17“^ 


-V 2 2^ - C, £ (1 - /„, )fc l2 + (a, - /„,a')v 2 + ft u 2 ^ 


2 ^ 


„„ k^\duv _ fc — <9v 
l) + 2Cc — v — — t-C$— mv — 
5 £ £ <?jy 




dy 

* , — <9t/ 

+C, (l-/wi)^i2 MV ^:- 


y,-/»ir’ + f-n‘' 2 


,£t/ 


+C 2 £(1 - 1^12 + ^12^22) “ 2 /wl 


uv 


k + 3v 2 / 2 


(C-37) 
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U^ + V^ = M{v + c e -v 2 )^ 

dx dy e Jdy 


~ C e i —uv 


dU 

dy 


~C £2 f + f w2 


r r e — du £ 2 ..££ 

L-UV-— + M N — 

y k dy k k 


(C-40) 
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C.5 Levy-Lees Transformation for Boundary-Layer Flows 

The body-oriented physical coordinates (x,y) and transformed coordinates (^,r\) are related 

by, 

and , (C-4.) 

where r n is the body radius, e stands for the freestream value and j=0, 1 for plane and 
axisymmetric flow, respectively. Equivalent to these expressions are the following: 


%(x) = ]p e u^r; , dx and 7](x,y) = ^^\ 


n 

( - \ 

f 

P 

J 

0 

\Pe) 


t’dy , 


(C-42) 


where t — r / r a is the transverse curvature. By using the chain-rule of calculus we can relate the 


derivatives of the two coordinate systems: 


d I 2j 

=Pe U eVe r o 


f .1 \ 


+ 


drj 


\dxj\dri 


f ^ \ 




(C-43a) 



. PAK 1 ' 

"p_^ 

(± \ 


1 

J-t 

\PeJ 



The dependent variables are transformed according to the following relations: 


F(£,rj) = -r- , 

u. 




V(^r?)= /f - 

Pe U 'Pe r „ 

7) = ^ , 


’dri} pvr J J J 


£(£,';) = ^lr 

U. 


(C-43b) 


(C-44) 


where the expression for V is extracted from the continuity equation assuming a stream-function 
exists. The equations governing the mean flow after the transformation are then given by, 

(C-45) 


■g- <9 F dV r n 
2 i H + Tn + F -° • 


2^ + v 3f * 




d£ dp dri 

For the k - e model the equations are transformed to, 


+ p(F 2 -Q-1) = 0 


(C-46) 
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2^F~ + V— - — 

dt, dp dr. 


t 2 ' L\ 


f n ^ 

J + tL 
\ °k j 


dK 


dl 7 


+ 2(3FK-t 2j Lp T 




(£ ~h £ ) 

+ ~ =q . 
p<pA j 


(C-47) 


2^F~ + V— - — 
< 9 £ dp dp 


t 2 ’ L\ 


( U ' 

/ + — 


<9£ 


<9 77 


+ 2(2[3 - 1)F£ 


£ f <9F V ^ /• 02 


~ C<#/# F ^ 


v<^y 


C f F 2 

+ -r~~ £ = 0 , 

p«p^' * 


where. 


£ = ^-e U ~C f 2 ^ P 
“ U 3 U 4 *’ Mt A ~1 


t-* 2 (2|) 


L(l + Q) 2 £ Uiptfu? 2 ' 

00 r e e r^e o 


(C-48) 


(C-49) 


The dimensionless parameters introduced are given by, 


£ - £ 

pJj„P,A 2l+l 

* u 

U e =— , 

U m 

O* ' — * » 

P.=&. 

(C-50a) 

r 

° A ’ 

- u 
Pe=~ , 

Pr 


(C-50b) 

W 

w- 

£- £ 



P„P r A 21 

P~P r A 2 ' 


(C-50c) 

The quantities. 




f-T 

T. 

L=P- P 

Pe P e 

and Q - 2< = du < 

*. d£ ’ 

(C-51) 


are the dimensionless temperature, density/viscosity and pressure gradient terms, respect.vely. 
The first two terms are introduced from the transformation and are zero and one, respectively, for 
isothermal, incompressible flows. The last term in (C-51) vanishes for flows with zero pressure 
gradient. The parameter 4 is a reference length and is set to 1 with no loss in generality. Similar 

transformations are also applied to the modeled transport equations of the near-wall Reynolds- 
stress models. 
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